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ABSTRACT 

Starting  with  single  server  queueing  systems,  we  find  a different 
way  to  estimate  the  diffusion  parameters.  The  boundary  condition  is 
handled  using  the  Feller's  elementary  return  process.  Extensive  comparisons 
by  asymptotic,  simulation  and  numerical  techniques  have  been  conducted  to 
establish  the  superiority  of  the  proposed  method  compared  with  conventional 
methods.  The  limitation  of  the  diffusion  approximation  is  also  investigated. 

When  the  coefficient  of  variation  of  interarrival  time  is  larger  than  one, 
the  mean  queue  length  may  vary  over  a wide  range  even  if  the  mean  and  variance 
of  interarrival  time  are  kept  unchanged.  The  diffusion  approximation  is 
applicable  under  the  condition  that  the  high  variation  of  interarrival  time 
is  due  to  a large  number  of  short  interarrival  times.  Case  studies  are 
conducted  on  2-stage  hyperexponential  distributions.  A similar  anomaly  is 
observed  in  two  server  closed  queueing  networks  when  the  service  time  of 
any  server  has  a large  coefficient  of  variation.  Again,  a similar  regularity 
condition  on  the  service  time  distribution  is  required  in  order  for  the 
diffusion  approximation  to  be  applicable.  For  general  queueing  networks, 
the  problems  become  more  complicated.  A simple  way  to  estimate  the  coefficient 
of  variation  of  interarrival  time  (when  the  network  is  decomposable)  is 
proposed.  Besides  the  anomalies  cited  before,  networks  under  certain  topologies, 
such  as  networks  with  feedback  loops,  especially  self  loops,  can  not  be 
decomposed  into  separate  single  servers  when  the  coefficient  of  variation 
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of  service  time  distributions  become  large,  even  if  the  large  variations 
are  due  to  a large  number  of  short  service  times.  Nevertheless,  the 
decomposability  of  a network  can  be  improved  by  replacing  each  server 
with  a self  loop  by  an  equivalent  server  without  a self  loop.  Finally, 
we  consider  the  service  center  with  a queue  dependent  service  rate  or 
arrival  rate.  Generalization  to  two  server  closed  queueing  networks 
where  each  server  may  have  a self  loop  is  also  considered. 
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1.  INTRODUCTION 


Recently,  considerable  effort  has  been  made  for  obtaining  approximate 
solutions  to  non-Markovian  queueing  models  using  the  diffusion  approxima- 
tion when  the  traffic  intensity  of  the  queueing  system  is  high.  The  ad- 
vantage of  diffusion  approximation  lies  in  the  fact  that  explicit  results 
can  be  obtained  for  relatively  complex  situations  where  the  only  possible 
alternatives  are  numerical  methods  or  simulation  experiments.  This 
greatly  extends  our  capability  in  modelling  practical  problems  . In  the 
past,  over  simplified  models  have  often  been  used  for  the  sake  of  mathe- 
matical tractability , and  the  predicted  ;-erformance  may  sometimes  be  quite 
different  from  the  actual  measured  performance. 

In  order  to  alleviate  the  difficulty  in^'clved  with  general  service 
time  distribution,  the  diffusion  approximation  replaces  the  discrete  .jump 
process  such  as  the  queue  size  process  by  a diffusion  process  which  is  a 
continuous  path  stochastic  process.  The  probability  distribution  of  the 
diffusion  process  which  satisfies  a partial  differential  equation  is  quite 
often  more  amenable  to  mathematical  analysis  than  that  of  the  .jump  process. 
However,  the  approximation  by  diffusion  process  requires  the  heavy  traffic 
assumption,  as  we  shall  see  in  Section  2. 

Based  on  central  limit  theorem,  Kingman  [8]  has  shown  in  his  treat- 
ment of  heavy  traffic  theory  that  the  waiting  time  distribution  is  as  an 
approximation  exponentially  distributed,  where  the  parameter  depends  only 
on  the  mean  and  variance  of  the  interarrival  time  and  service  time  dis- 
tribution, i.e.,  it  is  insensitive  to  the  detailed  form  of  the  distribu- 
tion, as  the  traffic  intensity  approaches  1.  The  diffusion  approximation 
based  on  the  same  idea  attempts  to  overcome  the  limitation  of  the  expo- 
nential model  by  considering  both  the  mean  and  variance  of  the  service 
time  and  Interarrival  time  distributions.  Newell  [11]  gives  an  extensive 
treatment  of  queues  with  time  dependent  arrival  rate  through  use  of  the 
diffusion  approximation  in  his  monograph.  Gaver  applies  the  diffusion 
approximation  method  to  waiting  time  in  a M/G/1  queue  [4].  Gaver  and 
Shedler  [2,3]  apply  this  technique  to  the  analysis  of  a multiprogrammed 
computer  system  modelled  as  a two  stage  cyclic  network.  Kobayashi  [10] 
considers  the  multi-dimensional  diffusion  approximation  as  a technique 
for  treating  general  queueing  networks.  Reiser  and  Kobayashi  [12] 
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study  the  accuracy  of  diffusion  approximation  techniques  and  propose  a 
way  to  treat  each  server  in  the  queueing  networks  separately.  Gelenbe 
[5,71  suggests  a different  way  to  handle  the  boundary  condition  of  the 
diffusion  process,  namely  using  the  Feller's  elementary  return  process 
[1] . In  [6],  Gelenbe  also  investigates  the  idea  of  decomposing  a queue- 
ing network  into  separate  single  servers.  An  application  of  the  diffu- 
sion approximation  to  analyze  the  performance  of  an  ALOHA-llke  system 
can  be  found  in  the  paper  by  Kobayashl,  Onozato,  Huynch  [17].  Kleinrock 
[9]  also  has  a tutorial  chapter  on  diffusion  approximation. 

Since  the  diffusion  approximation  to  single  server  queueing  systems 
serve  as  the  foundation  to  the  approximation  of  more  complicated  queueing 
networks,  we  will  start  with  single  server  systems,  then  advance  to  two 
server  closed  queueing  networks  where  each  server  may  have  a self  loop, 
and  finally  examine  the  problem  in  general  queueing  networks. 

In  Section  2,  we  propose  a new  way  to  estimate  the  diffusion  param- 
eters. Using  the  Feller's  elementary  return  process  [11  as  proposed  by 
Gelenbe  [5]  to  handle  the  boundary  condition,  the  approximate  mean  queue 
length  obtained  by  this  method  is  more  accurate  than  that  by  conventional 
methods  in  most  cases,  especially  when  the  coefficient  of  variation  of 
the  service  time  is  large.  In  Section  3,  we  analyze  the  as3miptotic  error 
in  mean  queue  length  by  our  method  and  two  other  widely  used  diffusion 
approximation  techniques  proposed  by  Kobayashl  [10]  and  Gelenbe  [5]  for 
the  M/G/1  and  E2/M/I  queueing  systems,  where  analytic  results  on  mean 
queue  length  are  available  in  closed  forms.  The  advantage  of  the  asymp- 
totic analysis  is  that  the  absolute  or  relative  errors  are  expressed  in 
terms  of  the  traffic  intensity  or  the  coefficients  of  variation  of  the 
service  time  and  interarrival  time  distributions.  This  provides  better 
insight  in  understanding  the  accuracy  of  various  approximation  techniques  . 
Kingman  [18]  has  found  a tight  upper  bound  for  the  mean  queue  length.  We 
also  analyze  this  upper  bound  for  reference.  It  is  interesting  to  see 
that,  in  the  M/G/1  system,  the  mean  queue  lengths  obtained  by  the  other 
two  methods  are  larger  than  the  Kingman's  upper  bound  when  the  service 
time  has  a large  coefficient  of  variation.  In  both  E^/M/1  and  M/G/1  sys- 
tems, our  method  yields  more  accurate  approximations.  In  fact,  in  the 
M/G/1  system,  the  mean  queue  length  obtained  by  our  method  is  exact,  and 
those  obtained  by  the  other  two  methods  have  an  error  term  on  the  order 
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of  C /2  or  (C  -l)/2,  \»4iere  C is  the  squared  coefficient  of  vari- 

s s s 

ation  of  the  service  time.  For  readers  who  are  not  familiar  with  asymp- 
totic analysis,  this  section  can  be  skipped  over.  In  Section  4,  simula- 
tions have  been  conducted  to  test  the  relative  p>erformance  of  our  method 
and  the  two  conventional  methods  for  more  general  queueing  systems  which 
include  the  E /E  /I  and  E /H_/l  systems.  Numerical  techniques  have  also 
been  employed  to  study  the  relative  performance  of  various  diffusion  ap- 
proximations for  the  Eg/M/l  and  D/M/1  systems.  Our  method  yields  more 
accurate  approximations,  except  in  the  E /E  /I  system.  In  the  E /E  /I 
system,  the  method  in  [4]  proposed  by  Kobayashi  has  better  performance 
than  ours.  These  comprehensive  and  systematic  comparisons  not  only  es- 
tablish the  robustness  of  the  proposed  method,  but  also  provide  valuable 
information  in  selecting  the  best  approximation  technique  for  the  spe- 
cific problem  at  hand.  In  Section  5,  we  use  the  system  to  illu- 

strate the  fact  that,  when  the  coefficient  of  variation  of  the  arrival 
process  is  larger  than  1,  the  mean  queue  length  may  vary  over  a wide 
range  even  if  the  mean  and  variance  of  the  interarrival  time  are  kept 
unchanged.  This  is  simply  because  the  coefficient  of  variation  of  the 
distribution  function  may  become  large  due  to  different  reasons.  We 
give  a reasonable  interpretation  to  this  phenomenon.  Since  two-stage 
hyperexponential  distribution  function  is  widely  used  in  computer  system 
modelling,  we  try  to  identify  the  range  of  the  parameters  of  the  hyper- 
exponential distribution  where  the  diffusion  approximation  can  be  ap- 
plied to  obtain  a fairly  accurate  estimation  of  the  mean  queue  length 
under  various  traffic  intensities.  The  data  included  in  that  section 
should  be  helpful  in  checking  the  applicability  of  diffusion  approxima- 
tion to  the  problem  at  hand.  The  H_/E  /I  and  H_/H_/l  systems  are  also 

<5  r ^ 2. 

examined.  In  those  cases  where  the  parameters  are  not  in  the  applicable 
range,  the  diffusion  approximation  may  be  used  to  estimate  a lower  bound 
of  the  system  performance. 

After  examining  the  single  server  queueing  system,  in  Section  6 we 
consider  a more  complicated  queueing  system,  the  closed  two  server  sys- 
tem which  is  often  used  to  model  the  computer  system  under  fixed  degree 
of  multiprogramming,  referred  to  as  the  CPU  and  DTU  model  [2,31.  The 
approximate  utilization  and  mean  queue  length  of  the  CF4J  are  very  close 
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to  the  simulation  or  exact  result  when  the  coefficients  of  variation  of 
the  service  time  distributions  are  small  [2,5].  When  this  condition 
does  not  hold,  the  diffusion  approximation  techniques  can  provide  a close 
approximation  to  mean  queue  length  and  utilization  only  under  restricted 
ranges  of  the  parameters  of  the  distribution  functions  and  in  other  cas- 
es, it  may  still  be  used  to  estimate  a lower  bound  of  the  performance  of 
the  system.  Then,  in  Section  7,  the  decomposition  problem  of  the  network 
of  queues  is  considered.  From  the  data  provided  in  Sections  3 and  4,  the 
accuracy  of  diffusion  approximation  to  a single  server  system  is  undoubt- 
edly very  good.  The  problem  on  decomposition  of  a queueing  network  into 
separate  single  server  systems  seems  to  be  how  to  estimate  the  coefficient 
of  variation  of  the  Interarrival  time  at  each  server,  so  we  can  take  the 
interactions  among  Interconnected  ser'^ers  into  account.  Two  different 
methods  have  been  proposed  to  estimate  the  coefficient  of  variation  of 
the  Interarrival  time  at  each  server  by  Reiser  and  Kobayashi  [12]  and 
Gelenbe  [6],  respectively.  Both  methods  lead  to  fairly  accurate  approx- 
imations. The  second  method  \dilch  tries  to  incorporate  the  effect  of 
idle  periods  on  the  coefficient  of  variation  seems  to  be  better  but  is 
more  complicated.  Here,  we  propose  a method  to  estimate  the  coefficient 
of  variation  of  the  Interarrlval  time  which  leads  to  similar  results  as 
the  second  method  by  taking  the  effect  of  idle  periods  into  consideration 
but  is  much  simpler  in  computation.  All  the  examples  given  by  the  pre- 
vious authors  to  demonstrate  the  accuracy  on  decomposing  queueing  net- 
works into  separate  single  servers  under  diffusion  approximations  are 
concentrated  on  the  situation  where  the  coefficients  of  variation  of 
service  time  and  external  interarrlval  time  distributions  are  not  large, 
mainly  less  than  or  equal  to  two.  Actually,  the  decomposition  technique 
is  not  always  feasible  when  the  coefficients  of  variation  of  the  service 
time  or  external  interarrlval  time  distributions  become  large.  An  exam- 
ple has  been  given  to  illustrate  this  anomaly  which  has  been  overlooked 
in  the  past.  Hence,  we  must  be  careful  on  the  decomposability  of  a queue- 
ing network.  Although  decomposability  is  an  inherent  property  of  the 
network  topology,  its  effect  magnifies  as  the  coefficients  of  variation 
of.  service  time  distributions  deviate  largely  from  one.  Nevertheless, 
decompKjsabil Ity  of  a network  can  be  Improved  by  replacing  each  server 
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with  a self  loop  by  an  equivalent  server  without  a self  loop.  Finally, 
in  Section  8,  we  consider  the  service  center  with  a queue  dependent  ser- 
vice rate  or  arrival  rate.  Generalization  to  the  closed  two  server 
queueing  network  is  also  considered. 
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2.  THE  DIFFUSION  APPROXIMATION  FOR  THE  G/G/1  QUEUE 


Consider  a single  server  queueing  system.  Let  t be  the  arrival 
th 

time  of  the  1 job  to  the  queueing  system  and  t’  be  the  departure  time 

of  the  1 job , where  0 < t.  < t < . • . , 0 < t ' < t ' < . . . , and  t I > 

X ^ X ^ X 

tj^,  i.e.,  the  queueing  discipline  of  the  system  is  first  come  first 
served  (FCFS).  Let  A(t)  and  D(t)  represent  the  cumulative  number  of 
arrivals  and  departures,  respectively,  up  to  time  t.  Denote  the  number 
of  jobs  in  the  queue  (including  the  job  in  service)  at  time  t by  Q(t), 
then 

Q(t)  = A(t)  - D(t) 


Assume  the  interarrival  time  U. 's  (=  t^  -t,  ,)  and  service  time  V. 's 

i i i-1  1 

are  independent  and  identically  distributed , respectively.  Furthermore, 
we  assume 


and 


Var  {U . } = cr^ 
1 a 


E(Vi)  = l 

Var {V  } = 

i s 


P 


A 

P 


where  p is  the  traffic  intensity  of  the  queueing  system. 

The  following  central  limit  theorem  for  renewal  processes  [14]  will 
be  used  in  later  discussion. 

Theorem.  If  T = (T  ) is  a renewal  process  (i.e.,  T -T  , are  inde- 
n ' n n-1 

pendent  and  identically  distributed)  for  which 


! 

I M s E(Tj^l  < 00 


6 


':T 


\ \ 


and 


Let 


= e|(Tj^  - < 00 

00 

N(t)  = A Irn  ) 

^ [0,t]  n 


wAiere 


I^(X)  = 
A 


if  X € A 


[o  if  X ^ A 


i.e.,  N(t)  is  the  number  of  T such  that  T < t.  Then 

n n — 


11,  P^N(t)  - t/M  < I ^ 


where 


t->oo  [ J 

^ loo  ^ 


du 


is  the  normal  integral. 

Now  let  us  introduce  the  definition  and  property  of  a diffusion 
process . 

Definitions . A diffusion  process  (x(t),t  >0}  is  a strong  Markov  pro- 
cess such  that 

(1)  B(x,t)  . u.  owt  ^ At)  - x(tUx(t)  . X) 

At^  0 

(2)  a(x,.)=  11.  oWt.iAO  - yt>/|x(t)  =x) 

At^O 

(3)  the  sample  path  is  continuous 

where  the  diffusion  parameters  p(X,t)  and  a(X,t)  are  called  the 
infinitesimal  mean  and  variance  coefficients,  respectively. 
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Let  P(Xp,X;t)  be  the  probability  density  function  of  the  diffusion 
process  X(t) , i.e., 

P(XQ,X;t)  dX  = p|x  < X(t)  < X + dx|X(0)  = Xq| 

then  P(XQ,X;t)  will  satisfy  the  following  differential  equation: 

^ P(Xo.X;t)  = i ^ a(x,t)  4 P(XQ,X;t)  - p(X,t)  PCX^.Xjt) 

Tidiich  is  called  the  Kolmogorov  diffusion  equation  or  Fokker-Plank  equa- 
tion [14]. 

Clearly,  as  t becomes  large,  the  renewal  counting  process  N(t) 
is  approaching  a diffusion  process  with 

«x,t)  = i 

and 

2 

a(x,t)  = -2- 

which  is  usually  referred  to  as  the  Wiener  process  with  drift. 

By  assumption,  the  arrival  process  is  a renewal  process,  hence  A(t) 

will  converge  to  a Wiener  process  with  Infinitesimal  mean  A and  infin- 

2 3 

Itesimal  variance  a A , as  t becomes  large.  The  problem  is  that  the 
interdeparture  time  is  not  independent  and  identically  distributed,  since 
the  Interdeparture  time  can  either  be  a service  time  or  the  sum  of  a ser- 
vice time  and  an  idle  period  of  the  server.  Hence,  the  departure  pro- 
cess is  not  a renewal  process.  But,  under  heavy  traffic  conditions,  i.e., 
as  p -> 1 , it  is  close  to  a renewal  process.  During  the  busy  period, 

D(t)  will  come  close  to  a Wiener  process  with  infinitesimal  mean  p and 

2 3 

variance  coefficient  as  t Increases,  provided  that  the  busy  pe- 

riod is  not  Interrupted.  Still  another  problem  is  that  A(t)  and  D(t) 
is  not  independent  since  Q(t)  = A(t)  -D(t)  > 0.  But,  when  Q(t)  is 
larger  than  zero,  we  have  a departure  process  independent  of  the  arrival 
process.  In  this  case,  Q(t)  behaves  like  a Wiener  process  \diich  is  a 
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diffusion  process  with  no  boundary  restriction  at  zero.  That  is  to  say, 
we  should  approximate  Q(t)  by  a diffusion  process  with  appropriate 
boundary  condition  at  zero  to  reflect  the  fact  that  Q(t)  can  never 
become  negative  and  there  is  an  idle  period  after  Q(t)  drops  to  zero. 

To  be  more  precise,  Q(t)  will  converge  to  a diffusion  process 
with  jjarameters 

p(X,t)  = 7\  - p 

2 3 2 3 

a(x,t)  = a A + a p p 

a s 


The  parameters  are  obtained  from  those  of  A(t)  and  D(t)  based  on  the 

fact  that  Q(t)  = A(t)  -D(t).  An  extra  factor  p appeared  in  the  second 

term  of  a(X,t)  is  used  to  reflect  the  fact  that  D(t)  has  a coefficient 
2 3 

of  variation  a^p  only  p of  the  time. 

Since  both  p(X,t)  and  a(X,t)  are  constant,  we  will  abbreviate 
them  as  g and  a,  respectively.  The  probability  density  function 
P(XQ,X;t)  of  Q(t)  will  satisfy  the  equation 

Si  = i ^ P(Xo.X,t)  - p ^ P(XQ,X,t) 

ox 

Let  P(X)  be  the  stationary  density  function  of  Q(t),  i.e., 

P(X)  d(X)  = p(x  < Q(t)  < X + dX}  as  t ^ oo 


For  the  stationary  case,  the  time  derivative  in  the  Fokker-Plank 
equation  is  set  to  zero.  So 


a 

2 


P(X)  = 0 


(2.1) 


Two  different  approaches  have  been  suggested  to  handle  the  boundary 
condition.  The  first  approach  is  to  treat  the  boundary  X = 0 as  a re- 
flecting boundary,  i.e.,  whenever  the  queue  becomes  empty,  it  is  reflected 
to  positive  immediately.  Though  the  queue  size  will  never  become  nega- 
tive, still  no  probability  mass  can  collect  at  X = 0.  Gaver  and  Shedler 
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[2,3]  and  Kobayashi  [10],  who  generalized  this  approach  to  queueing  net- 
work, have  managed  to  choose  the  appropriate  integration  constants  in  the 
solution  to  (2.1)  under  reflecting  boundary,  so  that  the  model  correctly 
predicts  the  stationary  probability  of  empty  queue.  The  second  approach 
proposed  by  Gelenbe  [5,6,7]  uses  Feller's  elementary  return  process  [1] 
Instead  of  the  diffusion  process  with  reflecting  boundary  to  approximate 
the  queueing  system.  This  is  a diffusion  process  with  boundary  to  which 
the  process  adheres  for  epochs  whenever  the  process  attains  a boundary; 
at  the  end  of  the  epoch  the  process  is  reinitialized  according  to  a fixed 
probability  density  function.  Gelenbe  [5]  first  solves  the  equation  when 
the  holding  time  on  the  boundary  has  exponential  distribution  and  later 
on  [7]  generalizes  it  to  any  probability  density  function  whose  Laplace- 
Stieltjes  transform  is  a rational  function  to  account  for  the  fact  that 
the  holding  time  in  the  boundary  in  general  is  not  exponentially  distrib- 
uted. Fortunately,  the  solution  under  the  general  distribution  depends 
only  on  the  first  moment  of  the  holding  time  distribution.  Thus,  the 
stationary  solution  is  identical  to  the  corresponding  solution  when  the 
holding  time  is  exponentially  distributed  with  the  same  mean. 

In  this  paper,  we  will  adopt  Gelenbe's  approach  to  handle  the  bound- 
ary condition  and  assume  exponential  holding  time  on  the  boundary,  since 
this  assumption  will  simplify  the  problem  and  lead  to  the  same  solution 
as  that  under  general  holding  time  distribution  [7].  The  advantage  of 
this  approach  is  that  it  can  be  extended  very  easily  to  handle  two  server 
closed  queueing  networks  or  finite  capacity  queues . Both  of  them  have  im- 
portant applications  in  computer  modelling.  The  only  problem  we  are  fac- 
ing is  that  the  mean  holding  time,  h,  at  the  boundary  X = 0 is  not 
known.  The  holding  time  at  the  boundary  X = 0 is , in  fact,  the  idle 
period  of  the  queueing  system.  From  queueing  theory  [23],  we  know  that 

h = E{idle  period) 

= A ^(1  - p)  E(n  ) 

where  E(n)  is  the  expected  number  of  jobs  being  served  in  each  busy 
p>erlod , and  furthermore 
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th 

Recall  V.  is  the  service  time  of  the  i customer  and  U.  is  the  in- 

1 th  th  ^ 

terarrival  time  between  the  i and  1-1  customer.  The  expression  for 

E(n)  can  be  simplified  to  l/(l-p)  when  the  interarrival  time,  U^,  has 

exponential  distribution,  l.e.,  the  system  is  M/G/1 . In  general,  it  can 

not  be  simplified.  We  will  use  the  conditions  that  the  integration  of 

probability  density  function  over  the  range  X > 0 should  equal  to  one 

to  obtain  an  estimation  of  the  holding  time,  h.  To  account  for  the  fact 

that,  after  an  arrival  to  the  empty  queue  occurs  the  number  of  customers 

in  the  queue  Jumps  instantaneously  to  one,  we  need  to  add  an  extra  term, 

-(1  - p)/h  6(X-1),  to  the  right  hand  side  of  (2.1)  and  an  extra  boundary 

equation  (2.3),  as  explained  below. 

Now  we  have  the  following  equations 

2 

I P(X)  - P Jj  P(X)  = - — B(X  - 1)  (2.2) 

c3x 

and 

“■  [f 

where  &(X  - 1)  is  a Dirac  density  function  concentrated  at  X = 1 and 
represents  the  probability  density  function  of  the  point  from  which  the 
diffusion  process  starts  once  again  immediately  after  a jump.  Notice 
(1  - p)/h  is  the  product  of  the  probability  at  the  boundary  and  the  rate 
of  jumping  back  from  the  boundary,  l.e.,  (1  - p)/h  represents  the  mean 

rate  of  jumping  back  to  (0,oo).  (o/2)(P/Px)  P(X)  - pP(X)  has  the  phys- 

ical interpretation  as  the  rate  of  flow  of  the  probability  mass  from  the 
region  (0,oo)  to  the  boundary  0.  This  explains  the  boundary  equation. 
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Similar  arguments  can  be  given  to  the  correction  term  in  (2.2).  For  a 
more  detailed  argument,  see  [5,7]. 

Let  us  denote 


where 


r 


-2(u  - A) 
2^3  2 3 

+ Ogii  p 


-2(1  - p) 

P(C  + C ) 
a s 


C 

a 


= a 


a 


C 

s 


2 2 


C and  C are  called  the  squared  coefficient  of  variations  of  the  in- 
a s 

terarrival  time  and  service  time,  respectively. 

Solving  the  differential  equation  (2)  with  boundary  condition  (3) 

and  lim  P(X)  =0,  we  get 
X-^  0 


P(X)  = <J 

/ (1  - p)  r.  -r-,  1 
(~~hp  [1  - e ] e 


To  compute  h,  we  use  the  fact  that 


(1  - P) 


+ 


P(X)  dX  = 1 


and  yield,  when  r<0  (i.e.,  p<l), 


As  we  can  see,  the  estimation  of  the  length  of  idle  period  by  diffusion 
approximation  is  only  exact  for  the  M/G/1  system. 

Finally,  we  face  the  problem  of  discretization  of  the  probability 
density  function  in  the  neighborhood  of  Integer  valued  points  X = i in 
oixler  to  approximate  , the  stationary  probability  of  finding  i cus- 
tomers in  the  queueing  system.  Usually  there  are  three  choices  for  n. : 
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3.  ACCURACY  ANALYSIS  OF  DIFFUSION  APPROXIMATION  BY  ASYMPTOTIC  TECHNIQUE 


In  this  section,  we  analyze  the  asymptotic  errors  of  the  mean  queue 
lengths  obtained  by  various  diffusion  approximation  techniques  for  the 
M/G/1  system  and  the  E^/iA/1  system  where  analytic  solution  of  the  mean 
queue  length  is  available  as  p-» 1 . Clearly,  it  is  an  important  require- 
ment for  an  approximation  technique  which  are  designed  to  handle  general 
distributions  for  service  time  or  interarrival  time  to  give  accurate  ap- 
proximation on  those  cases  for  which  solutions  are  known.  These  error 
analyses  certainly  have  important  implications  on  the  accuracy  of  the 
approximation  for  queueing  systems  whose  service  time  and  interarrival 
time  distributions  do  not  deviate  too  far  from  those  of  the  M/G/1  or 
E2/M/I  system.  The  advantage  of  asymptotic  analysis  is  that  we  can  ob- 
tain a closed  form  expression  for  the  error  term,  and  the  order  of  the 
error  can  be  clearly  expressed  in  terms  of  the  power  of  (1  - p)  which 
gives  us  a clear  picture  of  the  dependency  on  heavy  traffic  assumption. 
Various  diffusion  approximation  techniques  have  been  proposed  to  handle 
the  single  server  systems.  The  two  most  noteworthy  methods  are  proposed 
by  Reiser  and  Kobayashi  [10]  and  Gelenbe  [6],  respectively.  Since  there 
is  not  any  comprehensive  study  of  the  relative  accuracy  of  the  two  meth- 
(xis , we  will  analyze  not  only  the  proposed  method  but  also  the  two  meth- 
ods mentioned  above.  Our  method  improves  the  accuracy  in  both  cases. 

Tile  mean  queue  length  obtained  under  Kingman  heavy  traffic  approximation 
[H]  has  been  proved  to  provide  an  upper  bound  on  mean  queue  length  [18] . 
The  result  holds  for  0 < p < 1 and  improves  to  be  a tight  upper  bound 
as  p^l.  The  upper  bound  will  also  be  analyzed  for  comparison.  And  we 
shall  see  this  upper  bound  is  quite  tight  in  both  cases.  From  then  on, 
we  will  denote 

(1)  Method  P:  the  proposed  method 

(2)  Method  A:  the  diffusion  approximation  technique  proposed 

by  Kobayashi  [10] 

In  this  approximation  method,  we  have 

2 3 2 3 

a = CT  A + fT  ^ 
a s 
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P = n - A 


The  boundary  at  X = 0 is  treated  as  a reflection  boundary 
By  setting  the  probability  mass  at  the  origin  to  be  1 - p 
and  solving  the  Fokker-Plank  equation  (2.1),  the  following 
stationary  queue  length  distribution  and  mean  queue  length 
is  obtained. 

(a)  Stationary  queue  length  distribution 


rt 


0 


1 


P 


1 


p(l  - p) 


/\i_l 

P 


for  i > 1 


where 


(3.1) 


_ 2(l-p) 

C +pC 
^ s a 
p = e 


(b)  Stationary  mean  queue  length 


(3.2) 


(3)  Method  B:  the  diffusion  approximation  technique  proposed 
by  Gelenbe  [5] 

As  noted  earlier,  method  P follows  the  argument  in  method 
B to  handle  the  boundary  condition,  but  the  diffusion  pa- 
rameters in  the  two  methods  are  different.  The  diffusion 
parameters  in  method  B are 


2^3  2 3 

a = a A + a p 
a s 


/3  = M - A 


Hence,  the  stationary  queue  length  distribution  is 


«0  = 1 - c. 


p(pc  + c ) 

Si  s / r .,  . 

’'l  = 1(1  ~ p)  - 1 - r) 


(3.3) 


^ <=s>  lr„  -a,‘ 
”l  = -2(1~  P>'  ' ® ’ 


i > 2 


(3.3) 
Cont . 


_ 2g  ^ 2(u  - A) 

“a  2^3  2 3 

a A + a u. 
a s 


and  the  mean  queue  length  is 


/ pc  + C 


(3.4) 


Notice  the  Kingman's  upper  bound  on  mean  queue  length  [18]  is 


C + C p 

= 27— p'  ■ ♦ 0 


(3.5) 


As  we  shall  see,  in  the  M/G/1  system  the  mean  queue  length  obtained 
by  method  P is  exact  and  those  obtained  by  methods  A and  B have  an  abso- 
lute error  around  p(C^  - l)/2  and  pC^/2,  respectively.  Hence,  the 
performance  of  methods  A and  B will  degrade  as  the  coefficient  of  varia- 
tion of  the  service  time  Increases,  e.g.,  when  C equals  64  and  p 

s 

equals  0.8,  the  relative  errors  of  both  methods  are  around  25'^.  In  fact, 

under  heavy  traffic  condition  the  mean  queue  lengths  obtained  by  methods 

A and  B are  larger  than  the  Kigman ' s upper  bound  on  mean  queue  length 

when  C is  larger  than  3.  In  the  E_/M/l  system,  again,  the  mean  queue 
s ^ 

length  obtained  by  method  P is  more  accurate.  The  result  from  method  A 
is  very  close  to  that  from  method  P.  Examining  the  asymptotic  expres- 
sions for  the  mean  queue  lengths  of  both  methods,  we  find  their  differ- 
ence is  proportional  to  (1  - p) . The  mean  queue  length  obtained  by 
method  B is  less  accurate  and  is  quite  close  to  the  Kingman's  upper 
bound  on  mean  queue  length. 

We  will  first  prove  the  following  lemma,  which  is  the  foundation  of 
the  analysis  on  method  A. 
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Lemma  1.  The  asymptotic  mean  queue  length  obtained  under  method  A will 
satisfy 


E(Q^)  = 


P 


1 


A 

P 


where 


p(C  +pC  ) 
s a 


(1  - p) 


* 


6(C 


+ pc  ) 

a 


(1  - P) 


(3  .6 


P 


2(l-p) 

C +pC 
s a 


(3,7 


Proof . 


Using  the  Taylor  series  expansion  for 


X 

e , 


i ,e . , 


= 1 - X + It-  - It-  + o(x^ 


as  X 


we  obtain,  as  p^ 1 , from  (3.7) 


p = 1 - e 


2(l-p) 

C +pC 
s a 


= 1 


2(1  - p)  ^ 2(1  - p)^  ^ 4(1  - p)^ 

S (C  +pC  3(C  +pC 

s a s a 


2(1  -p) 

C + pc 
s a 


1 - p ^ 2(1  - p)^ 

S 3(C  -^PC 

s a 


(3.i 


Using  the  Taylor  series  expansion  of 


— = 1 + X 


2 3 

+ X + 0(X  ) 


1 


as  X ■>  0 


we  obtain,  as  p-> 1 , from  (3.8) 


p(C  +C  ) 
s a 

2 


(1  -P) 


+ 


P 

6(C  +PC  ) 
s a 


(1  - P) 


3.1  Mean  Queue  Length  for  M/G/1  System 

The  mean  queue  length  of  the  M/G/1  system  is  given  by  the  well  known 
Pollaczek-Khintchine  formula 


/ (1  + C )\ 

E(Q)  = p(l 

The  mean  queue  length  obtained  by  the  proposed  method  is  given  by 
(2.5),  Setting  C equal  to  one,  since  the  arrival  process  is  Poisson, 

£l 

we  find  that  the  mean  queue  length  given  in  (2.5)  is  exactly  the  Poliac- 
zek-Khintchine  formula.  That  is  to  say,  we  predict  the  mean  queue  length 
exactly  for  the  M/G/1  system.  Let  us  examine  the  asymptotic  performance 
of  method  A and  method  B in  this  case.  From  Lemma  1,  we  get  the  mean 
queue  length  under  method  A,  by  setting  = 1 , 

P(C  + p)  loo  / 2\ 

1 » - <»■  *1*  67c/—)-  " - P>  * »(<1  - P>  j 

and  the  absolute  error  is 
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E(Q^)  - E(Q)  = I <C^  - 1)  + 6(C  ♦ p)  « - P>  + 

s 


From  equation  (3.4),  we  get  the  mean  queue  length  under  method  B, 


by  setting  = 1, 


' *2(1 


+ c >■ 

s 

- p) 


and  the  absolute  error  is 


E{Qg]  - E[Q)  = — 


Since 


EfQg]  - E(QJ  = f (1  - P)  . 0 (1  - p)‘ 

s ' 


method  A and  method  B have  similar  performance  when  C is  large. 

s 

From  equation  (3.5),  we  get  the  Kingman's  upper  bound  on  mean  queue 
length  , by  setting  C = 1 , 


1 ^ ^ SP 

“2  (1  - p)  * ^ 


and  the  absolute  error  is 


E(Q^]  - E[Q]  = - (1  + p) 


As  we  can  see,  the  Kingman's  upper  bound  is  quite  accept ible  under 
the  M/G/1  case  and  differs  from  the  mean  queue  length  by  1/2(1  + p) . 
Although  the  absolute  errors  under  both  method  A and  method  B become 
large  as  the  coefficient  of  variation  of  the  service  time  increases,  the 
absolute  error  in  the  Kingman's  upper  bound  is  fixed.  For  3,  the 

mean  queue  length  obtained  by  method  A and  method  B is,  in  fact, 
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larger  than  the  Kingman's  upper  bound  when  the  traffic  is  high.  Let  us 

examine  the  relative  errors  in  both  methods  A and  B as  C becomes  large, 

s 

The  relative  error  for  method  A is 


E{Q^}-E(Q)  2 


- 1)  + 


e{q) 


6(C  + p) 

s 

2 (1  + C ) 

P s_ 

2 1 - p 


(1  - p)  + 


1-  0 ^(1  - p)^j 


y * p(l  + o) 


I!'  ■ (f  " 6(C^  " 


P (1  + C ) 

S 


0^(1  - p)^)^ 


1 - 


2(1  - 
p(l 


^ * o(<l  - p,^)) 


2(1  - p) 

P^(l  +C  ) 
s 


(c  - 1) 

s 


l6(C^  +p) 


C - 1 
s 

p(C  +1) 
s 


(1  -p)  + 


2^7p  - (5p  -6)C^  ~^^s) 

2 2 
p (1  + C ) (p  + C ) 
s s 


(1  -p)  + 


0 (d  - p)^j 


and  similarly  the  relative  error  for  method  B is 


E{Qg)  - E(Q) 


2C 


e(q) 


p(C  + 1) 

S 


(1  - p)  - 


p^d  + c 

s 


s , .2 

(1  - p)  + 


0^(1  - p)^j 


As  C becomes  large,  the  relative  errors  for  both  methods  A and  B 
s 

approach  d-p)/p.  Hence,  for  p equals  0.8,  the  relative  errors  are  25‘'i 
as  mentioned  earlier  (see  Table  3.1).  In  Tal les  3.1a  and  3.1b,  some  nu- 
merical comparisons  of  methods  P,  A,  and  B are  presented  for  = 0 , 

1/5,  1/4,  1/3,  1/2,  1,  2,  4,  8,  16,  32,  64,  and  128,  when  p = 0.9  and 
0.8,  respectively.  This  can  be  used  as  a check  for  the  correctness  of 
the  asymptotic  analysis. 
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Table  3.1a 


MEAN  QUEUE  LENGTH  FOR  M/G/1  SYSTEM  WHEN  p 


c 

s 

Method  P 
(Exact  Result) 

Method  B 

Method  A 

128 

523  .35 

580.95 

580 . 50 

64 

264.15 

292.95 

292.50 

32 

134.55 

148.95 

148.50 

16 

69.75 

76.95 

76 .50 

8 

37.35 

40.95 

40.50 

4 

21.15 

22.95 

22.50 

2 

13  .05 

13.95 

13.51 

1 

9.00 

9.45 

9.01 

1/2 

6.97 

7.20 

6.76 

1/3 

6.30 

6.45 

6.01 

1/4 

5.96 

6.07 

5.64 

1/5 

5.76 

5.85 

5.41 

0 

4.95 

4.95 

4.52 

Table  3.1b 

MEAN  QUEUE  LENGTH  FOR  M/G/1  SYSTEM  WHEN  p = 0.8 


c 

s 

Method  P 
(Exact  Result) 

— 
Method  B 

Method  A 

128 

207.20 

258.40 

258.00 

64 

104 . 80 

130.40 

130.00 

32 

53  .60 

66.40 

66.00 

16 

28.00 

34.40 

34  .00 

8 

15.20 

18.40 

18.00 

4 

8.80 

10.40 

10.00 

2 

5.60 

6 .40 

6.00 

1 

4.00 

4.40 

4.01 

1/2 

3 .20 

3 .40 

3 .02 

1/3 

2.93 

3 .07 

2.69 

1/4 

2.80 

2.90 

2.53 

1/5 

2.72 

2.80 

2.43 

0 

2.40 

2.40 

2.03 

3.2  Mean  Queue  Length  for  E„/M/l  System 

In  the  G/M/1  queueing  system,  the  stationary  distribution  of  number 
of  customers  in  the  system  is  given  by  [9] : 


Kq  = (1  - a) 


^ k-1 

= p(l  ~ a)  a k > 1 


(3.9) 


where  a is  the  unique  root  of 


a = A*(p  - |ia) 


(3.10) 


and  A*(s)  is  the  Laplace  Stieltjes  transform  of  the  interarrival  time 
distribution.  Furthermore,  the  mean  queue  length  is  given  by 


(3.11) 


When  the  interarrival  time  distribution  is  the  2-stage  Erlang  dis- 
tribution, we  can  get  a closed  form  solution  for  a,  i .e . , 


cr  - 4p  + 1 - V8p  + 1 


(3.12) 


Before  we  proceed  to  compare  the  asymptotic  behavior  of  various 
diffusion  approximation  techniques,  we  will  first  prove  the  following 


lemma . 


l£tm^.  The  steady  state  mean  queue  length  of  E /M/1  system  will  satisfy 


E{Q}  =1  p(l-p)"^  ^ (1-p)  + o^(l-p)2^ 


as  p -»  1 


From  (3.12), 


l-cr  = -(l-4p+  "Jsp  + 1) 
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Setting  Z = 1 - p,  after  simplification,  we  get 


1 - a = i («  - 3 +3J1  -|zj 
Using  the  Taylor’s  series  expansion  of  n/1  - X,  we  get 

1 - <,  = i (43  - 3 + 3(1  - I 2 - A 2"  - H.  Otz“>'jj 

Combining  (3.12)  and  (3.13)  together,  we  get 


(3.13) 


E{Q}  = 


Using  the  Taylor  series  expansion  of  (1-X)  , we  get 

E(Q)  - s (1  ♦ (I  z + H * *(h 


Finally,  substituting  1-p  for  Z,  we  get 


E{Q}  = I p(1  - p)“^  + ■*■  T§8  - P)  + 


0 ((1  - p)^^ 


From  equation  (2.5),  we  get  the  mean  queue  length  under  the  proposed 
method  P,  by  setting  C = 1/2,  C =1, 

3 S 

E{Qp}  = P + I P^(l  - P)~^ 
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BVom  equation  (3.4),  we  get  the  mean  queue  length  under  method  B, 

by  setting  C = 1/2,  C =1, 

& s 


= + I P^U  - p)  ^ + p 


and  the  absolute  error  is 


E{Qg}  - E[Q]  = I p - p(l  - p) 


Furthermore,  the  relative  error  is 


~ ® ^ 8 14  2 / 3) 

— TTqT 5 - P)  -gi  (1  - p)  Po((i  - B) 


0^(1  - B)^^ 


From  equation  (3.5),  we  get  the  Kingman  upper  bound  on  mean  queue 
length,  by  setting  C = 1/2,  C =1, 

Q S 


E{Q. 


+ P 


and  the  absolute  error  is 


E{a^)  - E(Q)  = i p)  - jIj  P(1  - 0)  + o((l  - 


From  the  above  analysis,  it  is  clear  that  the  mean  queue  length 
obtained  by  method  P has  minimum  absolute  error.  The  mean  queue  length 
obtained  by  method  A is  very  close  to  that  by  method  P,  since  the  dif- 
ference is  only  a first  order  term.  The  mean  queue  length  obtained  by 
method  B is  less  accurate.  In  fact , it  is  very  close  to  the  Kingman's  upper 

bound  on  the  mean  queue  length  under  heavy  traffic  condition.  The  rela- 

2 

tive  error  under  method  P will  approximately  be  2/9(1  - p)  - 7/81(1  - p)  , 
i.e.,  2.11,  for  p = 0.9  and  4.1')^,  for  p = 0.8,  as  can  be  checked  with 
Table  3.2.  Method  A has  similar  performance.  The  relative  error  of 
method  B will  be  approximately  8/9(1  - p)  - 14/81(1  -p)^,  i.e.,  8.7^ 
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for  p = 0.9  and  17.1^  for  p = 0.8,  as  can  be  checked  with  Table  3.2. 
Some  numerical  comparisons  of  methods  P,  A,  and  B are  presented  in  Table 
3.2  to  check  the  correctness  of  the  asymptotic  analysis  on  mean  queue 
lengths  and  their  errors  under  various  diffusion  techniques  for  the 
M/1  system. 
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p 


Exact  Result 


Method  P 


Method  B 


Method  A 


0.95 

14.331 

14.487 

14.962 

14.492 

0.90 

6.829 

6.974 

7.425 

6.985 

0.85 

4.327 

4 .463 

4.887 

4.477 

0.80 

3 .075 

3.200 

3.600 

3.219 

0.75 

2 .323 

2.438 

2,813 

2.460 

0.70 

1.820 

1.925 

2.275 

1.95 

4.  ACCURACY  ANALYSIS  OF  DIFFUSION  APPROXIMATION 
BY  SIMULATION  AND  NUMERICAL  TECHNIQUES 

In  the  previous  section,  we  analyze  the  accuracy  of  mean  queue 
lengths  obtained  by  various  diffusion  approximation  techniques  for  the 
M/G/1  system  and  the  E2/M/I  system  where  analytic  solution  of  mean  queue 
length  is  available.  Now  we  continue  the  accuracy  analysis  on  mean  queue 
lengths  obtained  by  various  diffusion  approximation  techniques  for  those 
queueing  systems  where  closed  form  expression  for  mean  queue  length  is 
not  available.  Either  simulation  or  numerical  technique  is  employed  to 
obtain  an  estimation  of  tne  mean  queue  length,  depending  upon  which  way 
is  more  convenient.  Since  simulation  is  only  a statistical  experiment 
and  its  convergence  is  slow  under  heavy  traffic  condition,  not  only  the 
point  estimation  but  also  the  95'/j  interval  estimations  are  included  to 
give  a better  feeling  on  the  accuracy  or  convergence  of  the  simulation. 

In  this  section,  we  will  concentrate  on  the  cases  where  the  coeffi- 
cient of  variation  of  the  arrival  process  is  less  than  1.  When  the  co- 
efficient of  variation  of  the  arrival  process  exceeds  1,  certain  anoma- 
lies might  happen,  as  we  shall  see  in  the  next  section.  We  will  choose 
Erlang  distribution  to  represent  the  Interarrival  arrival  time  distribu- 
tion since  its  coefficient  of  variation  is  less  than  one.  To  be  more 
specific,  the  squared  coefficient  of  variation  of  an  n stage  Erlang  dis- 
tribution is  equal  to  1/n  [9].  And  we  will  use  Erlang  distribution  and 
hyperexponential  distribution  to  represent  the  service  time  distribution 
with  coefficient  of  variation  less  than  and  greater  than  1,  respectively. 
A two  stage  hyperexp>onential  density  function  has  the  following  form: 

_ JL  _ JL 

~ M,  _ ^ " M„ 

W 1 (1  - U)  2 

— e + e 

Ml  M2 

where 

0 < w < 1 

Apparently,  it  is  the  combination  of  two  exponential  distributions  with 
mean  Mi  and  M2,  respectively.  The  probability  of  taking  the  first 
branch  is  w,  and  that  of  taking  the  second  branch  is  1 -(i'.  Let  us 
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■ 


r 

I 


iA 


assume  M,_,  < . After  simple  manipulation,  we  can  express  u'  and 

in  terms  of  , M,  and  C,  \^iiere  M and  C are  the  mean  and  squared 

coefficient  of  variation  of  the  distribution  function,  respectively. 


(M  - M^) 

- 2M2M  + (c  + 1) 


(4.1) 


M - (1  - w) 


(4.2) 


Furthermore , for  any  M and  C larger  than  1 , we  can  clioose  M._, 
arbitrarily  except  that  it  must  be  in  between  0 and  M.  The  latter  con- 
straint will  guarantee  that  the  w obtained  from  (4.1)  will  lie  in  be- 
tween 0 and  1,  and  the  obtained  from  (4.2)  will  be  positive.  Al- 
though various  combinations  of  w,  , and  will  lead  to  the  same 

mean  and  variance,  the  higher  moments  of  the  distributions  can  be  quite 
different.  Since  C can  be  chosen  arbitrarily,  we  can  get  any  value 
of  coefficient  of  variation  larger  than  one  by  using  two  stage  hyperex- 
ponential distribution  functions. 

As  jXJinted  out  earlier,  all  the  results  from  simulations  are  ex- 
pressed in  terms  of  95^j  confident  interval  estimations.  For  the  E /E  /^ 

n r 

system,  the  widths  of  the  confident  interval  are  less  than  4'  of  the 

point  estimations.  For  the  E /H„/l  system,  the  width  of  the  confidence 

Interval  grows  as  the  coefficient  of  variation  and  the  traffic  intensity 

increases.  Nevertheless,  even  in  the  worst  situation,  when  C =128  and 

s 

p = 0.80  in  Table  4.5  of  the  E„/H_/l  system,  the  interval  estimation 
of  mean  queue  length  is  203.7  ± 21  and  the  mean  queue  lengths  obtained 
by  methods  P,  B,  and  A are  206.4,  257.6,  and  257.2,  respectively.  Tlie 
superiority  of  methfxd  P is  apparent  in  this  case.  That  is  to  say,  in 
all  the  cases  where  simulations  are  used,  the  simulation  results  are  ac- 
curate enough  to  distinguish  the  relative  performance  of  various  diffu- 
sion approximations.  Otherwise,  numerical  technique  will  be  used. 

We  first  examine  the  E /E  /I  system.  Table  4.1  compares  the  results 

n r 

when  p = 0.85,  C = 0.5,  and  C = 1/2,  1/3,  1/4  and  1/5.  Method  A 
s a 


leads  to  the  most  accurate  result  within  3fo  relative  error.  Method  B is 

inferior  to  the  other  techniques  and  can  have  a relative  error  at  least 

up  to  17'^.  Method  P always  has  an  error  less  than  one  half  of  the  error 

in  method  B.  Table  4.2  contains  the  results  when  p=0.8,  C^=0.5, 

and  C = 1/2,  1/3,  1/4,  1/5.  Similar  patterns  are  again  observed.  As 
a 

we  shall  see,  the  case  where  both  coefficients  of  variation  of  interar- 
rival time  distribution  and  service  time  distribution  are  small,  i.e., 

the  E /E  /I  system  where  both  coefficients  of  variation  are  less  than  or 
nr 

equal  to  0 . 5 , is  the  only  case  where  method  P does  not  yield  the  best 
approxlmat ion . 

Let  us  now  examine  the  E /H_/l  system.  This  is  one  of  the  cases 

r 2 

where  method  P is  much  supreior  to  methods  A and  B.  Table  4.3  contains 

the  results  when  p=0.85,  C =0.5,  C =2,  4,  8,  16,  32,  and  64. 

s s 

The  relative  error  of  approximate  mean  queue  length  under  method  P is 
always  very  small  in  all  the  test  cases.  Methods  A and  B have  very  sim- 
ilar performance,  and  their  relative  errors  can  be  at  least  up  to  20 o 
when  p = 0.85.  If  we  look  at  the  table  more  carefully,  we  might  find 
that  the  absolute  errors  in  the  approximate  mean  queue  length  under  meth- 
ods A and  B are  very  close  to  (p/2)C^.  These  are  exactly  the  asymptotic 
absolute  errors  of  approximate  mean  queue  lengths  under  methods  A and  B 
in  the  M/g/1  system.  Again,  we  observe  the  robustness  of  method  P when 
the  coefficient  of  variation  of  service  time  distribution  is  large.  Ta- 
ble 4.4  compares  the  results  for  p = 0.85,  C = 1/3,  C =2,  4,  8,  16, 

3 S 

32,  and  64,  and  Tables  4.5  and  4.6  compare  the  results  for  C =0.5, 

3 

= 2,  4,  8,  16,  64,  and  128  when  p = 0.80  and  p = 0.75,  respec- 
tively. Similar  pattern  is  again  observed. 

Finally,  we  use  numerical  techniques  to  study  the  GI/M/1  system, 
namely  the  E^/M/1  system  and  the  D/M/1  system.  The  E,_,/M/l  system  has 
already  been  analyzed  in  Section  3,  and  we  pointed  out  in  (3.11)  that 
the  mean  queue  length  of  the  GI/M/1  system  is  p/(l  - n)  where  o is 
the  solution  of  the  equation  A*(u-pa)  = a.  For  the  E^/M/1  system, 
this  equation  is  a second  order  equation,  and  we  have  a closed  form  for 
the  root.  But,  for  the  E^/M/l  and  D/M/1  .ystems,  the  equations  arc 
third  order  and  transcendental  equations,  respectively;  so  we  have  to 
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Table  4.1 

MEAN  QUEUE  LENGTH  FOR  E /E  /I  SYSTEM  WHEN  p = 0.85 

n r 


c 

a 

c 

s 

Simulation 

Me  thod  P 

Method  B 

Method  A 

1/2 

1/2 

3.140  ± 0.039 

3.258 

3 .471 

3 .069 

1/3 

1/2 

2.694  ± 0.033 

2.857 

3 .069 

2.672 

1/4 

1/2 

2.478  + 0.024 

2.656 

2.869 

2.473 

1/5 

1/2 

2.347  ± 0.024 

2 . 536 

2.748 

2.355 

Table  4.2 

MEAN  QUEUE  LENGTH  FOR  E /E  /I  SYSTEM  WHEN  p = 0.80 

nr 


c 

a 

c 

s 

Simulation 

Method  P 

Method  B 

1— 

Method  A 

1/2 

1/2 

2.295  ± 0.022 

2.400 

2.600 

2.230 

1/3 

1/2 

1.986  ± 0.016 

2.133 

2 .333 

1.968 

1/4 

1/2 

1.840  ± 0.013 

2.000 

2.200 

1.838 

1/5 

1/2 

1.750  ± 0.013 
1 

1.920 

2.120 

1.760 
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Table  4.3 


i. 


MEAN  QUEUE  LENGTH  FOR  E2/H2/I  SYSTEM  WHEN  p = 0.85 


c 

s 

1 

M/M^ 

Simulation 

Method  P 

Method  B 

Method  A 

2 

2 

6.74  ± 0.19 

6.87 

7.72 

7.3 

4 

3 

11.52  ± 0.43 

11.69 

13  .39 

12.97 

8 

5 

20.77  ± 0.87 

21.32 

24.72 

24  .30 

16 

9 

40.63  ±2.54 

40.59 

47.39 

46.46 

32 

79.63  ± 6.46 

79.12 

92.72 

92.30 

64 

152.3  ± 14 

156.2 

183  .4 

183  .0 

1 


l 


Table  4.4 


MEAN  QUEUE  LENGTH  FOR  SYSTEM  WHEN  p = 0.85 


c 

s 

M/M^ 

Simulation 

Method  P 

Method  B 

Method  A 

2 

2 

6.26  ± 0.14 

6,47 

7.32 

6.90 

4 

3 

10.97  ± 0.35 

11.29 

12.99 

12.57 

8 

5 

20.25  ± 0,97 

20.92 

24.32 

23  .90 

16 

9 

40.20  ± 3.02 

40.19 

46.99 

46.56 

32 

17 

78,8  ± 6.9 

78,72 

92.32 

91.89 

64 

33 

152.0  ± 14 

155.8 

183  .0 

182  .6 

32 


Table  4.5 


MEAN  QUEUE  LENGTH  BOR  E2/H2/I  SYSTEM  WHEN  p = 0.80 


1 n 

i Simulation 

1 ^ 

Method  P 

Method  B 

Method  A 

4.67  ± 0.09 

4.80 

5.60 

5.21 

7.83  ± 0.22 

8.00 

9.60 

9.21 

14.11  ± 0.53 

14.40 

17.60 

17.20 

27.24  ± 1.39 

27.20 

33.60 

33  .20 

52.95  ± 3.02 

52.8 

65.6 

65.2 

102.4  ± 8 

104.0 

129.6 

129.2 

203 . 7 ± 21 

206.4 

257.6 

257.2 

Table  4.6 


MEAN  QUEUE  LENGTH  BOR  E2/H2/I  SYSTEM  WHEN  p = 0.75 


c 

s 

M/M^ 

Simulation 

Method  P 

Method  B 

Method  A 

■ 

2 

3.44  ± 0.05 

3.56 

4.31 

3 .95 

3 

5.67  + 0.12 

5.81 

7.31 

6.94 

8 

5 

10.08  ± 0.32 

10.31 

13  .31 

12.94 

16 

9 

19.27  ± 0.83 

19.31 

25.31 

24.94 

32 

17 

37.39  ± 1.92 

37.31 

49.31 

48.94 

64 

33 

73.02  ±4.73 

73.31 

97.31 

96.94 

128 

65 

146  ± 14 

145.3 

193.3 

192.9 

1 


3 2 2 3 

a - (2  +9p)  a + (1  +9p  +27p  ) a - 27p  = 0 for  the  E^/M/l  system 

and 

-(1-ct)u  „ , , 

e - cr  = 0 for  the  D/M/1  system 

In  Tables  4.7  and  4.8,  we  compare  the  mean  queue  lengths  obtained 
by  numerical  technique  with  those  by  various  diffusion  approximations 
for  the  Eg/M/1  and  D/M/1  systems,  respectively.  Again,  method  P is  the 
more  accurate  approximation  method.  The  mean  queue  length  obtained  by 
method  A is  very  close  to  that  by  method  P.  The  relative  error  in  meth- 
od B can  exceed  those  in  methods  A and  B by  25/,  in  some  cases.  Recall 
similar  phenomenon  appeared  in  our  analysis  on  the  Eg/M/l  system  in  Sec- 
tion 3. 


Table  4.7 

MEAN  QUEUE  LENGTH  FOR  E^/M/1  SYSTEM 


p 

Exact  Result 

Method  P 

1 

Method  B 

Method  A 

0.95 

12.775 

12.983 

13.458 

12.989 

0.90 

6.106 

6 .300 

6.750 

6.312 

0.85 

3.881 

4.061 

4.486 

4 .078 

0.80 

2.768 

2.933 

3 .333 

2.954 

0.75 

2.098 

2.250 

2.625 

2.275 

0.70 

1.650 

1.789 

2.139 

1.817 

Table  4,8 

MEAN  QUEUE  LENGTH  FOR  D/M/1  SYSTEM 


P 

Exact  Result 

Method  P 

Method  B 

Method  A 

0.95 

9.664 

9.975 

10.450 

9.983 

0.90 

4 .661 

4.950 

5.400 

4 .965 

0.85 

2.991 

3.258 

3.683 

3 .280 

0.80 

2.154 

2.400 

2.800 

2.427 

0.75 

1.651 

1.875 

2.250 

1.906 

0.70 

1.313 

1.517 

1.867 

1.551 

5.  ANOMALY  WHEN  THE  COEFFICIENT  OF  VARIATION  OF 
INTERARRIVAL  TIME  IS  LARGER  THAN  ONE 

In  this  section,  we  investigate  the  mean  queue  length  of  the  queue- 
ing system  where  the  coefficient  of  variation  of  the  interarrival  time 
is  larger  than  1.  We  will  use  the  H^/M/1  system  as  an  example  to  demon- 
strate the  anomaly  since  this  is  the  case  where  analytic  result  is  avail- 
able. In  the  M/G/1  system,  the  mean  queue  length  is  given  by  the  Pollac- 
zek-Khinchin  formula  (3.7),  and  it  depends  only  on  the  mean  and  variance 
of  the  service  time  distribution  and  the  mean  of  the  interarrival  time 
distribution.  This  seems  to  be  a support  of  the  robustness  of  the  dif- 
fusion approximation  which  only  utilize  the  means  and  variances  of  the 
interarrival  time  and  service  time  distributions  to  fit  the  parameters 
a and  3 of  the  Fokker-Plank  equation  and  neglects  the  effect  of  higher 
moments  of  the  distributions.  However,  after  examining  the  H2/M/I  sys- 
tem, we  will  find  that  higher  moments  of  the  interarrival  time  distribu- 
tion do  have  a drastic  effect  on  the  mean  queue  length  as  the  traffic 
intensity,  p,  deviates  from  1.  Since  the  two-stage  hyperexponential 
distribution  function  is  usually  used  \»hen  the  distribution  function  is 
required  to  have  high  coefficient  of  variation,  and  is  often  encountered 
in  computer  system  modelling,  we  will  further  analyze  the  regularity 
conditions  on  hyperexponential  distributions  for  the  diffusion  approxi- 
mation to  be  applicable.  That  is  to  say,  we  are  interested  in  identify- 
ing the  ranges  of  the  parameters  of  the  two-stage  hyperexponential  dis- 
tribution when  being  used  as  the  interarrival  time  distribution  such 
that  the  performance  of  the  queueing  system,  e.g.,  the  mean  queue  length, 
can  be  estimated  by  the  diffusion  approximation  accurately.  It  should 
not  be  too  surprising  that  the  ranges  will  shrink  as  p decreases  or 
Increases.  Recall  the  form  of  a two-stage  hyperexponential  distribution 
is  (u/M^)  + ((1-u)/M2)  , where  M2  < M^ . The  relations 

among  the  parameters  are  given  in  (4.1)  and  (4.2)  with  M = 1/X. 

Prom  Table  5. Id,  we  observe  that  by  choosing  different  w,  , and 
M2,  the  mean  queue  length  varies  over  a wide  range  from  624.62  to  23,18 
when  p = 0.95,  = 64.  Let  us  take  a closer  look  at  u,  , and  M,, 

at  two  extreme  cases: 


Case 

1; 

M^  = 32.818 

M^  = 0.01 

U = 3.018 

X 

10"^ 

Case 

2; 

M,  = 3151.0 

M„  = 0.99 

U = 3.175 

X 

io“® 

1 

2 

In  both  cases,  the  mean  of  interarrival  time  is  1,  and  the  variance  or 
the  square  coefficient  of  variation  is  64.  For  the  interarrival  time 
distribution  in  Case  2,  as  we  can  see,  M is  so  close  to  the  mean  in- 

/M  , 

terarrival  time,  the  first  exponential  density,  e 1,  has 

almost  no  effects  on  the  mean  interarrival  time  and  only  affects  the  va- 
riance. This  is  the  reason  why  the  mean  queue  length  under  the  interar- 
rival time  distribution  in  Case  2 is  almost  equal  to  the  mean  queue 
length  in  M/M/1  which  is  equal  to  19  when  p = 0.95.  That  is  to  say, 
although  having  C equal  to  64,  the  H_/M/l  system  in  Case  2 behaves 

Si  ^ 

very  much  like  an  M/M/1  system.  The  interarrival  time  is  generated  ac- 

/w  -X/M, 

cording  to  the  first  exponential  density  el  so  infre- 

quently that  we  may  ignore  it  when  evaluating  the  performance  of  the 
queueing  system.  But  in  Case  1,  AI2  is  so  close  to  zero,  the  first 
exponential  density  not  only  affects  the  variance  of  the  interarrival 
time  but  also  its  mean.  This  is  the  reason  why  the  mean  queue  length 
becomes  so  high,  624.62.  The  important  fact  to  realize  is  that  a large 
coefficient  of  variation  does  not  necessarily  mean  that  we  have  a lot 
of  short  interarrival  times,  and  the  fluctuation  is  quite  high,  as  in 
Case  1;  it  may  also  mean  there  is  few  very  long  interarrival  time  which 
makes  the  coefficient  of  variation  large  without  having  serious  effect 
on  the  fluctuation  of  the  system. 

Let  us  take  a look  at  the  mean  queue  length  predicted  by  the  diffu- 
sion approximations.  It  is  around  587  by  all  three  methods.  This  is 
quite  close  to  the  result  under  the  interarrival  time  distribution  in 
Case  1.  This  is  not  a surprise  since  mean  queue  length  calculated  by 
diffusion  approximation  is  fairly  close  to  the  Kingman's  upper  bound  on 
mean  queue  length  as  pointed  out  in  Section  3.  What  we  are  interested 
in  is,  can  we  decide  the  applicability  of  diffusion  approximation  to  the 
queueing  system  by  .lust  examining  the  parameters  of  the  distribution. 
Surely,  the  applicable  range  will  be  affected  by  p and  C . In  Tables 
5.1a  through  5. Id,  the  moan  queue  lengths  of  the  H^/U/1  systems  are  tab- 
ulated for  various  combinations  of  M^^/M,  M^/M,  w which  lead  to  the 
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Table  5.1a 

MEAN  QUEUE  LENGTH  WHEN  C 


1.5051 
1 . 5263 
1.5556 
1.6250 
1.7143 
1.8333 

2 .0000 
2.2500 
2.6667 

3.0000 
3 . 5000 

6.0000 

11.000 

50.998 


M /M 

U) 

Mean  Queue  Length 

1 

= 0.95 

p = 0.85 

0.01 

6.622  X lO"^ 

28.53 

8.50 

0.05 

6.435  X 10“^ 

28.48 

8.48 

0.1 

6.183  X lO"^ 

28.45 

8.45 

0.2 

5.614  X 10“^ 

28.39 

8.39 

0.3 

4.949  X 10“^ 

28.33 

8.34 

0.4 

4.186  X 10“^ 

28.26 

8.27 

0.5 

3.333  X 10“^ 

28.17 

8.19 

0.6 

2.424  X 10“^ 

28.06 

8.10 

0.7 

1.525  X 10~^ 

27.90 

7.96 

0.75 

1.111  X io“^ 

27.79 

7.87 

0.8 

7.407  X 10“^ 

27.63 

7.47 

0.9 

1.961  X 10“^ 

26.91 

7.29 

0.95 

4.975  X lO"^ 

25.76 

6.78 

0.99 

2.000  X lO"'* 

21.87 

5.97 

Table  5.1b 

MEAN  QUEUE  LENGTH  WHEN  C 


M^/M 

M^/M 

00 

Mean  Queu 

e Length 

p = 0,95 

p = 0.85 

4.5354 

0.01 

2.188  X lO"^ 

85.55 

25.47 

4.6842 

0.05 

2.050  X 10~^ 

85.35 

25.32 

4.8889 

0.1 

1.880  X 10~^ 

85.12 

25.12 

5.3750 

0.2 

1.546  X 10~^ 

84.67 

24  .67 

6.0000 

0.3 

1.228  X lo'^ 

84.11 

24  .12 

6.8333 

0.4 

-2 

9.326  X 10 

83,39 

23  .43 

8.0000 

0.5 

-2 

6.667  X 10 

82  .42 

22,51 

9.7500 

0.6 

4.372  X 10~^ 

81,02 

21,21 

12.667 

0.7 

2.507  X 10“^ 

78.79 

19.25 

15.000 

0.75 

1.754  X 10“^ 

77.03 

17.83 

18.500 

0,8 

1.130  X 10“^ 

74.49 

15.96 

36  .000 

0.9 

2.850  X lO"^ 

62.89 

10.53 

71 .000 

0.95 

7.138  X lO”'* 

45.94 

7.70 

351 .00 

0.99 

2.857  X lO"® 

23.16 

6 .00 

Table  5.1c 

MEAN  QUEUE  LENGTH  WHEN  C 


M^/M 

M^/M 

(jJ 

16.6565 

0.01 

5.947  X 10“^ 

17.3158 

0.05 

-2 

5.502  X 10 

18.2222 

0.1 

4.966  X 10“^ 

20.3750 

0.2 

3.965  X 10~^ 

23.1429 

0.3 

3.064  X lo"^ 

26.8333 

0.4 

2.270  X lO"^ 

32.0000 

0.5 

1.587  X lO"^ 

39.7500 

0.6 

-2 

1.022  X 10 

52  .6667 

0.7 

5.773  X lo"^ 

63 .0000 

0.75 

4.016  X lo'^ 

78 . 5000 

0.8 

-3 

2.574  X 10 

156.000 

0.9 

6.447  X lo"* 

311.000 

0.95 

1.613  X lo'^ 

1551.00 

0.99 

6.452  X lo'® 

Mean  Queue  Length 


p = 0 . 95  p = 0 . 85 


313  .65 
312.76 
311.81 
309.69 
306.98 
303 .47 
298.61 
291.18 
279.13 
269.58 
255.36 
187.90 
84  .62 
23  .43 


93, 

.37 

92 

.70 

91, 

.79 

(Ti 

00 

.68 

00 

.01 

83, 

.51 

00 

.65 

71 , 

.55 

60, 

.11 

51, 

.45 

39, 

.67 

Table  5. Id 

MEAN  QUEUE  LENGTH  WHEN  C =64 

a 


M^/M 

M^/M 

1 

w 

Mean  Queue  Length 

p = 0.95 

p = 0.85 

32.818 

0.01 

3.018  X 10“^ 

624  .62 

184.27 

34.158 

0.05 

2.785  X 10“^ 

616 .40 

182.61 

36  .000 

0.1 

2.507  X 10“^ 

614  .69 

180.73 

40.375 

0.2 

1.991  X 10“^ 

610.10 

176  .40 

46.000 

0.3 

1.532  X lo"^ 

604.62 

170.85 

53 . 500 

0.4 

1.113  X 10“^ 

596 . 79 

163 . 54 

64 .000 

0.5 

7.874  X 10“^ 

586.57 

153.35 

79.750 

0.6 

5.054  X 10“^ 

571.47 

138.28 

106 .00 

0.7 

2.850  X 10“^ 

545.90 

113 ,66 

127.00 

0.75 

1.980  X 10“^ 

525.69 

94  .63 

158.50 

0.8 

1.269  X 10~^ 

497.30 

67.99 

316 .00 

0.9 

3.174  X lO"'^ 

347.83 

14 . 10 

631.00 

0.95 

7.936  X 10“^ 

116.30 

8.02 

3151 .0 

0.99 

3.175  X 10~® 

23  .48 

6 .01 
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same  means  and  squared  coefficients  of  variation,  which  are  2,  8,  32, 
and  64,  respectively,  when  p = 0.95  and  0.85.  M^/M  and  U^/K  can  be 
viewed  as  the  normalized  means  of  the  first  exponential  and  second  expo- 
nential branches,  respectively.  As  we  can  observe  from  Table  5 .1 , as  ^2^^* 
increases,  the  mean  queue  length  decreases.  The  mean  queue  length  drops 
sharply  as  M^/M  approaches  1.  The  larger  the  coefficient  of  variation 
of  the  interarrival  time  is,  the  larger  the  variation  of  mean  queue 
length  can  be.  Nevertheless,  the  mean  queue  length  does  not  vary  too 
much  over  a substantial  range  of  the  value  of  M^/M.  This  gives  us  some 
hope  that  diffusion  approximation  may  be  applied  in  this  case.  In  Tables 
5.2a  and  5.2b,  the  mean  queue  lengths  of  the  H„/M/l  systems  when  C = 

2,  4,  8,  16,  32,  and  64  are  tabulated  for  p = 0.95  and  0.85,  respec- 
tively. The  exact  mean  queue  lengths  for  = 0.2,  0.5,  and  0.7  are 

also  included  in  Table  5.2a.  Similarly,  those  for  M2/M  = 0.1,  0.4,  and 
0.6  are  also  Included  in  Table  5.2b.  As  we  can  see,  various  diffusion 
approximations  lead  to  similar  results  and  they  provide  reasonable  ap- 
proximations to  the  analytic  results  under  those  M^/M  ratios.  Combin- 
ing the  results  from  Tables  5.1  and  5.2,  we  can  conclude  that,  for  p = 
0.95,  the  range  of  M2/M  where  diffusion  approximation  can  be  applied 
is 

'’2 

-M 

within  15'^  accuracy  for  < 64.  This  is  a very  conserved  bound.  When 
C is  close  to  1,  the  applicable  range  is  actually  larger  than  the  spec- 

3 

Ifled  range.  As  we  can  see  from  Table  5.1a,  for  the  case  C =2  and 

3 

p = 0.95,  even  when  M^/M  increases  to  0.95,  the  variation  of  mean  queue 
length  is  not  substantial,  and  the  mean  queue  length  from  diffusion  ap- 
proximation (see  Table  5.2a)  is  acceptable.  As  the  traffic  intensity, 
p,  decreases,  the  applicable  range  shrinks  very  quickly,  as  can  also  be 
observed  from  Table  5.1.  In  the  case  p = 0.85,  the  applicable  range 
is  around 

“2 

— < 0.6  to  0.65 
M — 

with  15‘i  accuracy  for  C <64. 

a — 
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16  154,4 

32  298.8 

587.6 


28.5 

46.5 

82.6 
154.8 
299.2 
588 


28,0 

46.1 

82.2 
154.4 
298,8 
587.6 


Mg/M  = 0,2 

Mg/M  =0.5 

28.4 

28.2 

47.2 

46.3 

85.1 

82.4 

160.7 

154.5 

309.7 

298.6 

610.1 

586,8 

Table  5.2b 

MEAN  QUEUE  LENGTH  FOR  H„/M/l  SYSTEM  WHEN  p = 0.85 


The  mean  queue  lengths  for  more  general  service  time  distributions 
under  input  are  not  available  analytically.  It  is  hard  to  do  a thor- 
ough Investigation  of  the  applicable  range  of  diffusion  approximation  in 
this  case  since  the  exact  mean  queue  length  can  only  be  estimated  through 
long  simulations  under  heavy  traffic  condition  or  tedious  numerical  tech- 
niques. Nevertheless,  we  select  certain  combinations  of  , w which 

fall  Inside  the  applicable  range  of  diffusion  approximation  to  ]I,_,/M/1 
systems  and  simulate  the  mean  queue  lengths  for  several  H^/H^/l  and  H^/ 

E /I  systems.  In  Table  5.3,  the  mean  queue  lengths  obtained  by  various 
r 

diffusion  approximations  are  compared  with  that  obtained  by  simulations 
for  the  H_/E  /I  system  when  p = 0,85  and  C =64,  32,  16,  8,  and  4. 

^ ^ 3. 

The  second  column  MgA  specifies  the  extra  degree  of  freedom  in  the 

Interarrival  time  distribution  and  also  is  the  criterion  we  use  to  test 

the  applicability  of  diffusion  approximation  in  the  H^/M/1  system.  The 

simulation  results  seem  to  be  quite  close  to  the  results  obtained  by 

various  diffusion  approximations.  All  three  diffusion  approximation 

techniques  yield  very  similar  results.  In  Table  5.4,  the  mean  queue 

lengths  obtained  by  various  diffusion  approximations  are  compared  with 

those  obtained  by  simulations  for  the  H^/H^/l  system  when  p = 0.85  under 

various  combinations  of  C and  C • The  interarrival  time  distribution 

a s 

is  assumed  to  have  the  same  form  as  before.  The  service  time  distribu- 
tion is  also  assumed  to  have  a similar  form  as  (u/M  ,)  e 4.  ((]  _ 

-X/M 

<^)/M  e ‘ with  mean  1/u.  The  second  column  is  the  same  as  that 
s2 

in  Table  5.3.  The  fourth  column  is  only  used  to  specify  the  extra  degree 
of  freedom  in  the  service  time  distribution.  The  results  under  v'arlous 
diffusion  approximations  are  again  quite  acceptable.  Furthermore,  meth- 
ods B and  A yield  very  similar  results,  and  the  result  from  method  P is 
somewhat  smaller.  Recalling  the  asymptotic  expressions  for  mean  queue 
lengths  in  Section  2,  we  know  that  the  mean  queue  length  obtained  by 
method  P is  smaller  than  those  obtained  by  methods  A and  B by  (p/2>  . 

This  is  Indeed  the  case  as  can  be  checked  from  Table  5.4.  Since  the  re- 
sults obtained  by  methods  A and  B are  closer  to  the  Kingman’s  upper 
bound,  the  mean  queue  length  obtained  by  method  P seems  to  cover  a wider 
range  of  the  parameters  of  the  distribution  within  accuracy.  Tliat 

is  to  say,  method  P seems  to  be  preferable  unless  M^A  is  very  close  to 
0 for  the  H2/H2/I  system. 
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^ ^ ^ 


Table  5.3 


MEAN  QUEUE  LENGTH  FOR  H^/K^/l  SYSTEM  WHEN  p = 0.85 


c 

a 

XM2 

Simulation 

Method  P 

Method  B 

Method  A 

64 

0.2 

163  ± 16 

156.2 

156.4 

156  .0 

32 

0.2 

83.5  ± 5.5 

79.1 

79.3 

78.9 

32 

0.5 

71.9  ± 6.3 

79.1 

79.3 

78.9 

16 

0.15 

43.5  ± 3.1 

40.6 

40.8 

40.4 

16 

0.5 

39.5  ± 2.6 

40.6 

40.8 

40.4 

8 

0.3 

22.6  ± 1.1 

21.3 

21.5 

21.1 

8 

0.45 

21.5  ± 0.8 

21.3 

21.5 

21.1 

4 

0.3 

12.2  ± 0.4 

11.7 

11.9 

11.5 

4 

0.5 

11.5  ± 0.4 

11.7 

11.9 

11.5 

Table  5.4 


MEAN  QUEUE  LENGTH  FOR  H^/H^/l  SYSTEM  WHEN  p = 0.85 


c 

a 

AM.^ 

c 

s 

f^’s2 

Simulation 

Method  P 

Method  B 

Method  A 

64 

0.2 

64 

0.2 

312  ± 33 

309 

336 

336 

0.2 

32 

0.5 

237  ± 19 

232 

246 

24  5 

32 

0.2 

64 

0.6 

230  ± 20 

232 

259 

259 

16 

0.3 

16 

0.3 

74.8  ± 5.0 

77.9 

84.7 

84  .3 

8 

0.75 

4 

0.75 

25.0  ± 1.4 

29.8 

31.5 

31  .0 

B 

0.6 

4 

0.6 

18.4  + 1.0 

20.1 

21.8 

21  .4 

1 

0.75 

4 

0.75 

18.5  ± 0.6 

20.1 

21.8 

21  .4 

From  then  on,  we  will  use  the  term  "tjTje  a"  hyperexponential  dis- 
tribution to  denote  the  hsrperexponential  distribution  having  the  property 
that  M^/M  is  not  "close"  to  1.  We  will  be  vague  on  the  exact  range  of 
the  parameters  that  tirpe  A hyperexponential  distributions  must  satisfy, 
l.e.,  we  only  define  it  in  a qualitative  way,  not  in  a quantitative  way. 
For  more  complicated  queueing  systems,  the  exact  range  of  the  parameters 
of  type  A hyperexponential  distributions  such  that  the  diffusion  approx- 
imation can  be  applied  to  obtain  a reasonable  estimate  of  mean  queue 
length,  utilization,  etc.  may  be  different,  dep>ending  on  the  network 
topology  and  traffic  Intensity  at  the  server.  But,  as  long  as 
not  far  from  zero,  the  approximation  should  be  applicable.  Furthermore, 
we  will  use  the  term  "type  b"  hyperexponential  distribution  to  represent 
the  other  extreme,  l.e.,  M2/M  is  close  to  1. 

Although  hyperexponential  distributions  have  been  used  throughout 
the  section  as  interarrival  time  distributions,  we  expect  the  results 
should  hold  for  other  general  distributions.  That  is  to  say,  as  long  as 
the  high  variation  of  interarrival  time  is  due  to  a large  number  of  short 
Interarrlval  times  instead  of  a few  very  long  Interarrlval  times,  diffu- 
sion approximation  should  be  applicable.  Similarly,  the  terms  type  A and 
type  B distributions  are  used  to  denote  the  two  types  of  distributions 
with  high  coefficients  of  variation,  respectively. 


6.  CLOSED  TWO  SERVER  SYSTEMS  (CPU/DTU  MODEL) 


After  completing  the  discussion  on  diffusion  approximation  to  single 
sei*ver  systems,  we  now  consider  applying  diffusion  approximation  to  closed 
two  server  systems,  as  shown  in  Fig.  6.1.  In  computer  system  modelling, 
the  closed  two  server  system  is  often  used  to  represent  the  CPU  (central 
processing  unit)  and  DTU  (data  transfer  unit)  operating  under  fixed  de- 
gree of  multiprogramming.  From  then  on,  we  will  call  the  two  servers  CPU 
and  DTU,  respectively.  This  model  has  been  analyzed  by  Gaver  and  Shedler 
[21  using  diffusion  approximation  with  reflecting  boundaries  and  usual 
way  to  estimate  a and  p when  the  CPU  service  time  which  is  the  time 
between  page  faults  under  this  particular  interpretation  has  exponential 
distribution,  i.e.,  when  its  coefficient  of  variation  is  equal  to  1. 
Ijater  on,  Gaver  and  Shedler  [3]  use  Wald's  identity  to  fit  the  ratio  of 
29,/U  and  analyze  the  same  system  when  the  CPU  service  time  has  hyper- 
exponential distribution,  i.e.,  when  its  coefficient  of  variation  is 
larger  than  1.  Gelenbe  [51  also  analyzed  this  system  using  the  Feller's 
elementary  return  process  and  the  usual  way  to  estimate  diffusion  param- 
eters . 


Fig.  6.1.  CLOSED  TWO  SERVER  SYSTEM  (CPU-DTU  SYSTEM). 


The  anomalies  in  Section  5 lead  us  to  suspect  the  same  kind  of  prob- 
lems may  exist  in  the  two  server  closed  queueing  networks  if  one  or  more 
servers  have  hyperexponential  service  time  distributions.  This  is  simply 
because  the  departure  process  of  any  server  is  the  arrival  process  of  the 
other  server.  Before  we  proceed  to  investigate  the  anomalies,  let  us 
first  examine  the  various  diffusion  approximation  techniques  on  the  two 
server  closed  queueing  network  lust  mentioned  in  some  detail.  We  will 
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concentrate  on  the  CPU  system  and  try  to  estimate  the  utilization  or  the 
mean  queue  length  of  the  CPU  system.  This  measure  is  of  practical  im- 
portance since  it  can  be  used  as  an  indicator  of  the  stationary  system 
performance.  Assume  the  fixed  degree  of  multiprogramming  is  M.  Let 
P(X)  be  the  stationary  probability  density  function  for  the  diffusion 
process  in  (0,M),  where  X is  the  number  of  programs  in  the  CPU  queue 
and  M-X  is  the  number  of  programs  in  the  DTU  queue.  Furthermore,  let 

us  assume  the  mean  service  time  of  the  CPU  and  DTU  are  l/p  and  l/>\, 

2 2 

and  the  variance  of  the  CPU  and  DTU  are  d and  a , respectively. 

s a 

The  first  diffusion  approximation  technique  which  we  are  going  to 
examine  is  the  method  proposed  by  Gaver  and  Shedler  [2].  In  this  method, 
we  have 

2, 3 2 3 

a = cr  A + a u 

a s 


P = H - A 


as  usual.  Let  F(X)  be  the  stationary  distribution  function  of  the 
queue  length  at  CPU.  Imposing  a reflection  boundary  at  X = 0 and  nor- 
malizing the  probability  mass  between  0 and  M to  one,  we  get 


F(X)  = 


1 - A e 


rX 


1 - A e 


rM 


where 


2p  _ 2(^  - A) 

a ~ 2^3  2 3 

CT  A + CT  p 

a s 


by  solving  the  Fokker-Plank  equation  (2.1).  The  unknown  constant  A is 
chosen  such  that 


lim  F(0)  = 1 - p 

M CXI 


After  simple  manipulation,  we  get 


A = p 
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Hence , 


F(X) 


P e 
P e 


rX 

rM 


We  will  refer  to  this  method  as  method  G1  later  on. 

The  second  diffusion  approximation  technique  is  again  due  to  Gaver 
and  Shedler  [3]  to  handle  the  case  -where  the  coefficient  of  variation  of 
the  CPU  service  time  is  larger  than  one.  Let  G(S)  be  the  Laplace- 
Stleltjes  transform  of  the  CPU  service  time  distribution  and  H(S)  be 
the  Laplace-Stielt jes  transform  of  the  DTU  service  time  distribution. 
Furthermore,  let  S*  be  the  positive  solution  of 

G(-S)  H(S)  = 1 


By  Wald's  identity,  we  get 

r = ^ = In  G(S*) 


Following  the  same  argument  as  the  previous  method,  we  get 


F(X) 


1 


1 


A e 
A e 


rX 

rM 


Now,  using  the  fact  that  the  long  run  input  rate  to  the  CPU,  AF(M- 
1),  must  equal  to  the  long  run  output  rate  from  the  CPU,  p(l  -F(0)), 
we  get 


, -r(M-l)  -rM 

1 + p e - e 

We  will  refer  to  this  method  as  method  G later  on. 

Finally,  we  consider  the  method  proposed  In  Gelenbe  [51  with  a 
slightly  different  argument  from  [5]  making  it  coherent  with  Section  2. 
Again,  we  refer  to  this  method  as  method  B. 
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Now  there  are  two  boundaries:  0 and  M.  When  the  number  of  pro- 

grams in  the  CPU  queue  is  M,  the  DTU  is  idle  and  no  new  arrival  will 
occur.  It  is  assumed  that  the  CPU  queue  length  will  .jump  to  M-1  after 

an  exponential  holding  time  with  mean  h„.  This  is  similar  to  the  way 

M 

the  boundary  at  X = 0 is  handled  in  Section  2.  The  boundary  at  X = 0 
is  still  handled  in  the  same  way  as  in  Section  2,  where  the  mean  holding 
time  is  assumed  to  be  h^ . As  noted  earlier,  restricting  the  holding 
time  distribution  on  the  boundaries  to  be  exponentially  distributed  does 
not  mean  to  impose  any  restriction  on  the  service  time  distribution  of 
the  CPU  or  DTU.  It  is  .just  a conceptual  help  for  us  to  handle  the  bound- 
ary conditions.  Now  we  are  facing  a more  serious  problem  than  we  en- 
countered in  the  GI/G/1  system.  Here,  not  only  the  two  mean  holding 
times  at  the  boundary  X = 0 and  X = M are  unknown,  but  also  the 
probabilities  of  having  empty  queue  and  full  queue  are  unknown.  As  we 
shall  see,  the  boundary  equations  can  only  be  used  to  solve  two  unknowns, 
i.e.,  we  must  find  approximate  values  for  two  of  the  four  unknown  param- 
eters. Since  the  probability  of  empty  queue  is  directly  related  to  the 
CPU  utilization,  a quantity  of  ma,jor  concern,  we  will  try  to  find  reas- 
onable estimations  for  h.,  and  h , the  mean  holding  times  at  the  bound- 

1 M 

arles,  and  hope  that  the  errors  in  these  estimations  will  have  minor  in- 
fluence on  other  quantities  of  interest.  Recall  the  GI/G/1  system  in 
Section  2 where  the  approximate  holding  time  at  X = 0 obtained  by  the 
diffusion  approximation  is  equal  to  1/X,  which  is  the  mean  interarrival 
time.  Hence,  we  will  set  h^  equal  to  l/>. , the  mean  service  time  of 
the  DTU  since,  when  the  DTU  is  busy,  the  mean  interarrival  cime  to  the 
CPU  is  equal  to  the  mean  service  time  of  the  DTU.  By  similar  argument, 

we  will  set  h„  equal  to  1/u,  the  mean  service  time  of  the  CPU. 

M 

At  steady  state,  we  have  the  following  equations: 


f 


,! 


fl 


1 

I iJ 

I i; 


/l  ^ 

' M 


)X 


P(X)  - 3 P(x)  = 


-AM^6(X  - 1)  - - M + 1) 


lim  ^ a ~ P(X)  - 3P(X) 
I X ^ 0 ^ 


= AM, 


lim  I a P(X)  - 3P(X) 
X M 


t| 


] 
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where  &(.•)  is  the  Dirac  density  function  and  6(X  - 1)  and  6(X-M+1) 
represent  the  probability  density  function  of  the  point  from  which  the 
diffusion  process  starts  once  again  immediately  after  a jump  from  the 
boundaries  0 and  M,  respectively;  and  are  the  probability 

masses  concentrated  on  the  lower  boundary  and  upper  boundary,  respec- 
tively. 

Furthermore,  we  have  the  following  boundary  conditions: 


lim  P(X)  = lim  P(X)  = 0 
X^  0 X^  M 


Solving  the  above  equations,  we  get 


AM, 


1 -rX. 

- (1  - e ) 


0 < X < 1 


-r  rX 

P(X)  = ( (e  - 1)  e 1 < X < M - 1 

P ~ “ 


(6.1) 


Am, 

\“3 


2 ^^r(X-M)  M-1<X<M 


where 


and 


^1  r(M-l) 

”2  = — ® 


r 

a 


Also,  using 


/■ 

•'n 


P(X)  dX  + M^  + M^  = 1 


we  get 


M^  = (1  - p)  1 


2/  r(M-l)\“^ 
P (e  j 


(6.2) 
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A more  detailed  derivation  of  the  result  is  given  in  Gelenbe  [5], 
We  will  further  consider  the  problem  of  discretization  of  the  probabil- 
ity density  function  in  the  neighborhood  of  integer  valued  point  x = i 
in  order  to  approximate  , the  stationary  probability  of  having  i 
jobs  in  the  system.  The  following  way  of  discretization  is  proposed. 


2 < i < M - 2 


After  simplification,  we  get 


In  the  G/G/1  system,  we  find  that  the  accuracy  of  diffusion  approx- 
imation can  be  improved  by  defining  the  diffusion  parameter  3 as  + 

pgC^  where  g is  set  to  p,  the  traffic  intensity  or  the  utilization 
of  the  server.  For  the  closed  two  server  system,  the  utilization  of  each 
server  is  not  apparent.  From  experimental  results,  it  seems  to  be  that 
setting  g equal  to  A/p  for  > 1 and  1 otherwise  may  improve 
the  accuracy.  We  will  denote  this  method  as  method  P. 

After  examining  the  various  diffusion  approximation  methods,  now 
let  us  consider  the  case  where  the  coefficient  of  variation  of  the  ser- 
vice time  distribution  at  the  CPU  is  large,  as  is  often  the  case.  If  the 
service  time  distribution  at  the  DTU  is  exponentially  distributed,  the 
system  is  analytically  tractable  since  it  can  be  viewed  as  a M/G/1  queue- 
ing system  with  finite  waiting  room  [19].  We  summarize  the  result  for 
stationary  mean  queue  length  as  follows  . 

M 

1 < i < M 

i M i — 


and 


where 

is  the  stationary  probability  that  the  queue  length  is  equal  to 
i when  the  capacity  of  the  waiting  room  is  M 

is  the  stationary  queue  length  distribution  of  the  M/G/1  system 

Hence,  we  will  assume  the  service  time  at  DTU  has  exponential  distribu- 
tion and  use  the  analytic  result  to  analyze  the  accuracy  and  applicabil- 
ity of  diffusion  approximation.  The  answer  to  the  following  three  ques- 
tions are  of  major  Interest:  (1)  Does  the  mean  queue  length  or  the 

utilization  of  the  CPU  vary  if  type  B hyperexponential  distribution  in- 
stead of  type  A hyperexponential  distribution  is  used  for  the  service 


M , ^ 

’^M  = ^ ^ 


/ M-1 

p(l  - ^ Jt 


1 - 


^ 1 
i=0 
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time  distribution  of  CPU,  (2)  Does  the  diffusion  approximation  give  a 
reasonable  approximation  to  the  mean  queue  length  and  utilization  when 
type  A hyperexponential  distribution  is  used  for  service  time  distribu- 
tion of  CPU,  (3)  Does  the  service  time  of  the  CPU  often  have  type  A 
hyperexponential  distribution,  i.e.,  we  are  more  interested  in  the  ap- 
plicability to  computer  system  modelling. 

The  hyperexponential  distribution  of  the  CPU  service  time  is  as- 
sumed to  have  to  form  (w/M^^)  e ^"^^^1  + ((1-(jj)/M2>  e ''^'^*2  with  mean  1/u 
where  M.^  < . From  Tables  6.1  and  6.2,  we  see  that  the  mean  queue 

length  and  utilization  change  as  the  parameters  of  the  hyperexponential 
distribution  change.  Similar  phenomenon  on  CPU  utilization  in  the  M'G'' 
1/N  system  is  observed  by  Price  [24],  As  M^u  decrease,  i.e.,  the  num- 
ber of  requests  having  short  CPU  interval  increases,  the  analytic  re- 
sults become  very  close  to  the  results  obtained  by  both  diffusion  ap- 
proximation methods. 

Again,  we  see  the  diffusion  approximation  gives  a good  approximation 
of  the  performance  under  type  A hyperexponential  service  time  distribu- 
tion. We  also  expect  the  results  can  be  generalized  to  more  general  dis- 
tributions. That  is  to  say,  diffusion  approximation  will  be  applicable 
if  the  large  variation  of  service  time  is  due  to  a lot  of  short  service 
times.  If  it  is  due  to  a few  long  service  times,  diffusion  approximation 
can  only  be  used  to  obtain  a lower  bound  of  the  performance. 

The  third  question  is  hard  to  answer  in  general.  In  [31,  three 
sets  of  data  on  the  CPU  utilization  and  the  mean  and  variance  of  the 
CPU  service  time,  which  is  the  time  between  page  faults,  are  given.  The 
mean  and  variance  of  the  CPU  service  time  are  gathered  from  actual  pro- 
gram data.  The  results  on  CPU  utilization  are  obtained  by  trace  driven 
simulations  of  that  queueing  system.  The  DTU  service  time  is  assumed 
to  be  constant  to  account  for  the  average  access  time  along  with  the 
time  to  transfer  a page  of  information.  In  Tables  6.3,  6.4,  and  6.5, 
we  compare  the  results  obtained  by  three  different  diffusion  approxima- 
tion techniques,  i.e.,  methods  P,  G,  and  B,  with  the  analytic  result 
obtained  by  semi-Markov  analysis  [221.  The  parameter  M2 ' ® 
perexponent ia 1 distributions  taken  in  [31  are  all  less  than  four  tenths 
of  their  means,  respectively,  so  they  should  be  type  A hyperexponential 
d istribut Ions  . 
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IZATION  OF  CPU  WHERE  CPU  HAS  SERVICE  TIME  DISTRIBUTION  WITH 


Table  6.3 


l/[l  = 17026, 


CPU  UTILIZATION 
10 


a = 0.39780  X 10 
s 


= 3682,  1/?  = 20,000 


Number  of  Jobs 
in  the  System 


2 

3 

4 

5 

6 

7 

8 


Semi-Markov 

Method  P 

Method  G 

Meth<xl  B 

0,5316 

0.4934 

0.5993 

0 ,4887 

0 .5548 

0.5223 

0.6667 

0,5141 

0.5752 

0.5475 

0.7064 

0.5367 

0.5935 

0 . 5696 

0,7326 

0.5568 

0,6098 

0 . 5892 

0.7511 

0.5749 

0.6245 

0.6067 

0,7650 

0,5912 

0.6379 

0.6223 

0.7757 

0.6060 

0.6500 

0.6364 

0.7842 

0.6195 

0,6611 

0.6491 

0.7911 

0.6318 

Table  6.4 
CPU  UTILIZATION 
^2 


2 q 

l/\i  = 4871,  cr  = 0,26492  X 10'  , M = 1929,  lA  = 20,000 


Number  of  Jobs 
in  the  System 

Semi-Markov 

Method  P 

Method  G 

Method  ) 

2 

0.2216 

0.2169 

0.2216 

0.2022 

3 

0.2286 

0.2285 

0.2313 

0.2077 

4 

0.2333 

0.2350 

0,2361 

0.2124 

5 

0.2366 

0.2387 

0.2388 

0.2165 

6 

0.2388 

0.2408 

0.2404 

0.2201 

7 

0.2403 

0,2420 

0.2415 

0.2231 

8 

0 .2413 

0,2426 

0.2422 

0.2258 

9 

0,2420 

0.2430 

0.2426 

0.2281 

10 

0.2425 

0.2432 

0.2429 

0.2301 
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Table  6.5 

CPU  UTILIZATION 
2 10 

10735.  a = 0.12313  X 10  . M = 2953.  lA  = 20.000 


Number  of  Jobs 
in  the  System 


Semi-Markov  Method  P Method  G Method  R 


0.3704 
0.3887 
0.4045 
0 .4183 
0.4304 
0.4410 
0.4505 
0.4588 
0.4663 


0.4076 

0.3803 

0.4249 

0.4281 

0.4147 

0.4579 

0.4449 

0.4368 

0.4764 

0.4587 

0 .4544 

0.4882 

0.4702 

0 .4685 

0.4964 

0 .4798 

0.4799 

0.5024 

0.4879 

0.4893 

0.5070 

0 .4948 

0 .4970 

0.5106 

0 . 5006 

0.5033 

0.5136 

As  expected,  the  diffusion  approximations  are  very  close  to  the 
analytic  result,  and  method  P seems  to  have  better  overall  performance 
among  the  three  diffusion  approximations.  In  Tables  6.6a,  6.6b,  and 
6.6c,  we  compare  all  of  them  with  the  results  obtained  by  trace  driven 
simulations  in  [3],  and  the  results  are  again  very  close.  Hence,  we  can 
say,  at  least  in  that  computer  environment,  the  assumption  of  having  type 
A hyperexponential  service  time  distribution  is  very  reasonable.  Lewis 
and  Shedler  [20],  based  on  a statistical  analysis  of  actual  computer 
program  address  traces,  presented  a semi-Markov  model  for  the  point  pro- 
cess of  page  exceptions.  Although  the  detailed  stochastic  structure  of 
that  model  is  more  complicated,  it  does  have  the  property  that  it  con- 
sists of  a large  number  of  short  interfault  time  similar  to  that  of  type 
A hyperexponential  distribution. 

For  the  case  where  the  coefficient  of  variation  of  CPU  is  small 
(<1),  methods  P and  B are  equivalent.  In  [5],  comparison  of  methods  B 
and  G1 , with  results  obtained  by  semi-Markov  analysis,  shows  that  both 
methods  are  very  accurate,  and  method  B is  a little  bit  better. 

The  model  in  Fig.  6.1  can  also  be  generalized  to  include  a self  loop 
at  each  server,  as  shown  in  Fig.  6.2.  To  account  for  the  effect  of  the 
self  loop,  we  treat  each  server,  including  its  self  loop,  as  a single 
entity  and  consider  the  effect  of  the  self  loop  as  an  internal  interac- 
tion which  is  transparent  to  other  parts  of  the  system.  That  is  to  say, 

* 

we  will  replace  the  server  with  a self  loop  by  an  equivalent  one  without 

a self  loop.  The  interdeparture  time  of  the  equivalent  server  is,  in 

fact,  the  service  completion  time  seen  by  DTU . The  mean  and  variance  of 

2 

the  service  time  of  the  equivalent  CPU  are  (1  -0)/u  and  (cr  /(l-f>)  + 

2 2 ® 

<9/(u  (1-0)  )),  respectively.  The  two  quantities  are  derived  below. 


Fig.  6.2.  CPU-DTU  MODEL  WITH  SELF  LOOPS. 
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* in  the  sense  that  queue  length  distribution  and  departure  rate  to  the 
other  server  are  preserved. 


[ '' 


y 


Let  N be  the  number  of  service  completions  at  CPU  in  between  jobs 
arriving  at  the  DTU  (including  the  last  one).  Then,  N is  a geometri- 
cally distributed  random  variable  with  mean  l/(l-0)  and  variance  9/(1  - 
2 

9)  . Let  X be  the  random  variable  which  represents  the  service  time 
of  CPU  and  Y be  the  random  variable  which  represents  the  service  time 
of  the  equivalent  CPU  without  self  loop.  By  assumption, 

1 2 

e(x]  = — and  Var[x)  = a 
p.  s 

Clearly , 

E{Y)  = E(e[y|n}] 

= E{NE(X)} 

1 

u(l  - 9) 

Using  the  Identity  for  conditional  variance  [141,  we  get 


Var{Y)  = ErVar(Y|N)]  + Var[E(Y|N)] 
= E[N  Var{x)]  + Var{NE{x}] 

= E^Ndf}  + Var{^| 

2 


Similarly,  we  can  derive  the  mean  and  variance  of  the  service  time  of 

equivalent  DTU  without  self  loop.  The  mean  and  variance  of  the  service 

2 2 2 

time  are  equal  to  (l-ij/)/X  and  (ct  /(1-j/)  + ii/(u  (1  - \1/)  )),  respec- 

a 

tively.  After  obtaining  the  means  and  variances  of  equivalent  servers 
without  self  loop,  we  reduce  the  model  to  the  original  closed  two-server 
queueing  model . 

When  both  stages  have  exponential  service  times,  the  model  in  Fig. 
6.2  is  equivalent  to  that  in  Fig,  6.1  with  service  rates  p(l  -9)  and 
>,(1  - ,')  at  CPU  and  DTU,  respectively.  The  forms  of  the  service  time 
distributions  do  not  change  in  this  case. 
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7.  GENERAL  QUEUEING  NETWORKS 


Finally,  let  us  consider  applying  diffusion  approximation  to  ana- 
lyze the  performance  of  queueing  networks.  Kobayashi  [10]  proposed  that 
queueing  processes  of  a general  queueing  network  be  approximated  by  a 
vector-valued  diffusion  process.  The  interactions  among  different  que- 
ueing processes  are  explicitly  considered  in  the  diffusion  equations  in 
terms  of  the  variance-covariance  matrix.  The  Joint  queue  length  distri- 
bution is  expressed  in  a product  form  of  the  marginal  queue  size  distri- 
butions. This  solution  form  suggests  us  to  treat  each  queue  separately 
by  properly  taking  into  account  the  interaction  among  different  queues 
[12] . From  our  analytic  and  experimental  data  presented  in  the  previous 
section,  we  know  the  accuracy  of  diffusion  approximation  is  extremely 
good  on  single  server  system  except  for  certain  pathological  cases  given 
in  Section  5.  So,  the  success  on  decomposing  queueing  network  into  sep- 
arate single  server  systems  solely  relies  on  whether  we  can  find  a good 
estimation  of  the  coefficient  of  variation  of  the  interarrival  time  dis- 
tribution or  the  coefficient  of  variation  of  the  interdeparture  time 
distribution  at  each  server  such  that  the  correlations  among  servers  are 
not  ignored  after  decomposition.  Two  different  methods  to  estimate  the 
coefficient  of  variation  of  the  interdeparture  time  at  each  server  have 
been  proposed  by  Reiser  and  Kobayashi  [12]  and  Gelenbe  [6],  respectively. 
The  method  proposed  by  Gelenbe  tries  to  take  into  account  the  effect  of 
idle  period  on  the  coefficient  of  variation  of  interdeparture  time  \diich 
is  neglected  by  the  other  method  and  seems  to  lead  to  better  results . 
Nevertheless,  this  method  is  more  complicated  in  the  sense  that  matrix 
inversion  is  involved.  In  this  section,  we  propose  a simpler  way  to  es- 
timate the  coefficient  of  variation  of  the  interdeparture  time  distribu- 
tion than  that  by  Gelenbe,  yet  the  effect  of  idle  period  is  taken  into 
account.  The  values  obtained  by  both  methods  are  close  to  each  other. 
Furthermore,  all  the  demonstrating  examples  given  by  the  previous  auth- 
ors to  show  the  accuracy  on  decomposing  queueing  network  into  separate 
queues  using  diffusion  approximation  are  under  the  condition  that  the 
coefficient  of  variation  is  not  large,  mainly  loss  than  or  equal  to  2. 

Wo  present  an  anomaly  which  has  been  overlooked  in  the  past,  i.o.,  cer- 
tain network  topology  can  only  be  decomposed  into  separate  single  servers 
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when  all  the  service  times  and  external  interarrival  times  are  "nearly" 
exponentially  distributed.  When  the  coefficients  of  variation  of  some 
of  the  service  time  or  external  interarrival  time  distributions  deviates 
further  from  1 or  the  traffic  intensity  decreases,  the  decomposition  of 
this  kind  of  queueing  network  will  not  be  feasible  if  we  still  adopt  the 
conventional  way  to  estimate  the  coefficients  of  variation  of  interde- 
parture times  or  diffusion  parameters. 

Let  us  first  examine  the  method  in  [121  proposed  by  Reiser  and  Ko- 
bayashi  to  estimate  the  coefficient  of  variation  of  the  interarrival 
time.  In  their  treatment,  the  coefficient  of  variation  of  the  service 
time  is  taken  to  be  the  coefficient  of  variation  of  the  interdeparture 
time  as  a simple  approximation.  Furthermore,  the  departure  processes 


from  different  servers  are  treated  as  independent  renewal  processes. 


. th 


After  considering  the  fact  that  the  departure  process  of  the  i server 


are  only  active  p.  percent  of  the  time,  where  i . is  the  utilization 

th  ^ ^ i 

of  the  i server,  they  obtain  the  following  expression  for  C , the 

^ th 

squared  coefficient  of  variation  of  the  interarrival  time  of  the  i 


server ; 


n 

^ ^ / [(C.  - 1)  P..  + l1  A.P.. 

“ in  L .1  ,n  J .1  .li 


i ,i=o 


) P. . +1 


(7.1) 


Furthermore , C . 


.1 


the  squared  coefficient  of  variation  of  the  interde- 


th 


parture  time  of  the  ,j  server,  is  approximated  by 


C . = C 


(7.2) 


where 


.th 


is  the  arrival  rate  to  the  ,i  server  for  i 1 


is  the  external  arrival  rate 


is  the  squared  coefficient  of  variation  of  tlio  service  time 
distribution  in  the  server 


ii 


is  the  routing  probability  that,  after  departing  from  tlie  i 
server,  the  ,iob  will  i<jin  the  i^*'  server  queue,  for  ,j  ■ 1 


til 


Oi 


is  the  probability  that  the  external  arrival  will  loin  the 


server  queue 


’>3 
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Notice 


equations 


{\] 


is  the  solution  of  the  following  system  of  linear 


Oi 


j=i 


A.P., 
J Ji 


(7.3) 


A detailed  interpretation  of  this  equation  can  be  found  in  [13].  For 
open  queueing  networks,  (7.3)  provides  unique  solution  to  ] . In 

this  section,  all  the  queueing  networks  considered  are  assumed  to  be 
open  queueing  networks.  Extension  to  closed  queueing  networks  follows 
the  treatment  in  [12]. 

Gelenbe  [6]  argued  that  Reiser  and  Kobayashi  [12]  did  not  put  into 

consideration  of  the  effect  of  idle  period  on  the  variance  of  the  inter- 

th 

departure  time  and  assumed  that  the  interdeparture  time  of  the  i ser- 
ver, is  a service  time  with  probability  or  an  interar- 

rival time,  , plus  a service  time  with  probability  (1-p^). 

By  straightforward  manipulation,  we  get 


= e|s^J  + (1  - p^)  ('Kl  + 2e|a^| 


for  1 < i < n 


and  by  definition 


e|t^|  = A~^(l  + C^) 


for  1 < 1 < n 


Combining  the  two  equations  together,  we  get 


Ci^i 


for  1 < i < n (7.4) 


Finally,  making  the  assumption  that  the  number  of  arrivals  forms  a 
renewal  process,  we  get 
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^ F(c.  - 1)  P..  + i]  A.p  = - (a7^)  ) 

L .1  -u  J j .n  i y I i)  \ ^ / / 


for  1 < i < n 


After  simplification,  we  get  the  following  system  of  linear  equations: 


C.  = p^(C^  + 1)  + (1  - p.)(2p.  + 1 + A"’- 
1 1 s 1 \ 1 . 


^ [(C.  - 1)  P.. 


+ 1 


A.P..  - 1 

J u/ 


for  1 < i < n 


(7.5) 


The  assumption  that  is  equal  to  service  time,  S^,  with  pi-ob- 

ability  p.  or  an  interarrival  time  plus  a service  time,  A^ 

probability  (1-p.)  is  exact  only  for  Markovian  queueing  netvvoi’ks  . In 

^2  2 
Markovian  queueing  networks,  E[a.]  is  equal  to  2/A..  Hence,  a sim- 

^ 2 ^2 

pier  approach  is  to  further  approximate  e{A^]  by  2/7>^;  then,  we  get 
from  (7.1)  that 


That  is  to  say,  we  can  now  directly  express  in  a closed  form  expres- 

sion, and  the  necessity  of  solving  a system  of  linear  equations  has  been 
eliminated.  The  value  obtained  by  both  methods  are  very  close,  as  we 
shall  see. 

Alternatively,  we  may  consider  the  problem  in  the  following  way. 
Let  (z  be  the  set  of  service  centers  whose  customers  after  service  com- 
pletion may  go  to  service  center  A for  further  service.  Since  the  ar- 
rival rate  to  each  service  center  in  the  open  queueing  network  is  known, 
we  first  unfold  the  network  but  retain  the  connections  from  those  service 
centers  in  f to  A.  If  A is  in  f,  i .e . , there  is  self  loop  at  A, 
a duplication  of  A is  used  to  replace  the  self  loop.  Tlien,  wc  apply  a 
Poisson  input  to  each  service  center  in  t'  with  the  same  arrival  rate 
as  before  and  then  calculate  its  coefficient  of  variation  of  ^nterrlopai'- 
ture  time.  Finally,  treating  each  departure  process  as  an  inde;.'-  ..lent 
renewal  process,  we  can  obtain  the  coefficient  of  variations  of  the  in- 
torarrival  time  at  service  center  A in  this  case.  Tlie  coefficient  of 
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variation  obtained  will  be  used  as  an  approximation  of  that  of  the  orig- 
inal network  and  is  exactly  the  same  as  substituting  (7.6)  into  (7.1). 
For  example,  in  the  network  given  in  Fig.  7.1,  the  appropriate  subnet- 
work for  evaluating  the  coefficient  of  variation  of  the  interarrival 
time  of  the  service  centers  1 and  2 are  given  in  Figs.  7.2a  and  7.2b, 
respectively . 


Fig.  7.1.  OPEN  TWO  SERVER  QUEUEING  MODEL. 

We  now  use  the  queueing  network  shown  in  Fig.  7.1  to  compare  the 
' s , the  squared  coefficients  of  variation  of  the  interdeparture  time 
at  each  server,  by  our  straightforward  method,  Gelenbe's  method,  and 
Reiser  and  Kobayashi's  method.  Both  authors  [6,12]  have  used  this  net- 
work to  demonstrate  the  accuracy  of  their  approximations.  In  Table  7.1, 
the  column  under  method  P*  contains  the  results  of  the  proposed  method, 
and  the  column  under  method  B*  contains  the  results  of  the  method  pro- 
posed by  Gelenbe  [6],  and  the  column  under  method  A*  contains  the  results 
of  the  method  proposed  by  Reiser  and  Kobayashi  [12].  As  we  can  see,  all 
three  methods  yield  the  same  answer  \dien  the  network  is  a Markovian  que- 
ueing network.  The  estimates  of  the  coefficient  of  variation  of  the  inter- 
departure time  obtained  by  our  method  and  Gelenbe's  method  are  always  very 
close,  as  it  should  be  since  both  methods  try  to  incorporate  the  effect 
of  idle  period  on  the  coefficient  of  variation  of  the  interdeparture 
time.  However,  our  method  is  much  simpler  in  computation. 
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(a)  First  server 


Xe„ 


(b) 


Second  server 


Fig.  7.2.  APPROXIMATE  NETWORK  CONFIGURATION  FOR 
ESTIMATING  THE  COEFFICIENTS  OF  VARIATION  OF  IN- 
TERARRIVAL TIMES  FOR  THE  NETWORK  IN  FIG.  7.1. 


t 


Table  7.1 


' THE  SQUARED  COEFFICIENTS  OF  VARIATIONS  OF  INTERDEPARTURE  TIMES  (C^ , 

Cg)  OF  THE  QUEUEING  NETWORK  IN  FIG.  7.1  WITH  = 0.9,  = 0.84 

3> 

:4~ 

(a)  = 9^  = 0.5 


V 


i , 

I 


j' 


Let  us  apply  the  diffusion  approximation  techniques  to  analyze  the 
computer  communication  network  in  Fig.  7.3.  The  network  has  the  same 
topology  as  the  communication  network,  CIGALE,  within  CYCLADES  [28]  which 
is  a general  purpose  computer  network  being  Installed  in  France.  All  the 
terrestrial  links  are  assumed  to  be  full  duplex.  The  numbers  on  the 
terrestrial  links  represent  servers  and  their  queues.  Thus,  3 refers  to 
the  server  which  transfers  messages  from  node  C to  node  A and  2 refers 
to  the  server  which  transfers  messages  in  the  opposite  direction.  Traf- 
fic moving  in  the  two  opposite  directions  along  the  same  link  is  assumed 
to  be  non  interfering . Each  station  receives  external  traffic  which  forms 
a Poisson  process.  We  also  assume  that  each  message  arriving  from  out- 
side to  station  i has  equal  probabilities  of  having  any  of  the  other  1 
stations  as  its  final  destination.  The  routing  algorithm  of  the  networks 
is  assumed  to  be  fixed  and  will  be  described  later.  All  the  above  as- 
sumptions about  the  terrestrial  network  have  been  adopted  by  Gelenbe  [6] 
in  modeling  the  CYGALE  network  under  packet  switching. 


Fig.  7.3.  COMPUTER  COMMUNICATION  NETWORK. 
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Let  cH^  be  the  channel  capacity,  the  number  of  packets  that  can  be 
transmitted  per  second,  of  link  1.  The  channel  capacity  of  each  link 
is  indicated  in  Table  7.2.  The  fixed  routing  algorithm  is  summarized  in 
Table  7.3.  The  routes  which  are  not  shown  in  Table  7.3  are  the  links 
which  directly  connect  the  source  stations  and  destination  stations.  The 
number  of  packets  contained  in  each  message  is  assumed  to  be  geometri- 
cally distributed  with  mean  five.  The  external  arrival  rate  at  each 
node  is  tabulated  in  Table  7.4. 

The  performances  of  the  network  under  both  message  switching  and 
packet  switching  are  analyzed.  In  Table  7.5a,  mean  queue  lengths  at 
each  server  under  message  switching  by  various  approximation  methods  and 
simulation  are  tabulated.  Again,  our  method  is  denoted  by  method  P, 
Gelenbe's  method  is  denoted  by  method  B,  and  Reiser  and  Kobayashi's 
method  is  denoted  by  method  A.  The  simulation  results  are  presented 
with  95'^  confidence  Intervals.  In  Table  7.5b,  the  corresponding  squared 
coefficients  of  variation  of  interdeparture  time  at  each  server  by  vari- 
ous  methods  are  tabulated.  Both  methods  P and  B lead  to  similar  re- 
sults on  the  squared  coefficients  of  variation  of  interdeparture  times. 
The  difference  in  mean  queue  by  methods  P and  B in  Table  7.5a  is  mainly 
due  to  different  ways  being  employed  in  estimating  diffusion  parameters. 
The  minor  difference  in  estimating  squared  coefficient  of  variation  of 
interdeparture  time  has  very  little  effect.  As  we  can  see,  method  P 
leads  to  better  approximation.  In  Table  7.6a,  the  mean  queue  lengths 
under  packet  switching  obtained  by  methods  P,  B,  A and  simulation  are 
tabulated.  In  Table  7.6b,  the  corresponding  squared  coefficient  of  va- 
riation of  interdeparture  time  at  each  server  by  various  methods  are 
tabulated.  Not  only  the  squared  coefficients  of  variation  of  interde- 
parture  time  but  also  the  mean  queue  lengths  obtained  by  methcxi  P and 
method  B are  very  close  to  each  other.  Furthermore,  the  mean  queue 
lengths  obtained  by  both  methods  are  very  close  to  the  simulation  re- 
sult. As  we  can  see  from  Tables  7.5a  and  7.6a,  method  P provides  very 
accurate  approximations  in  both  cases,  where  other  methods  can  provide 
very  accurate  approximations  in  only  one  of  the  two  cases. 

Finally,  we  consider  the  decomposabil Ity  problem  of  general  queue- 
ing networks.  From  our  previous  analyses  in  Sections  5 and  6,  we  expect 
that,  if  the  service  time  distributions  of  some  of  the  intermediate 
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Table  7.2 

CHANNEL  CAPACITY  OF  EACH  LINK 
IN  THE  TERRESTRIAL  NETWORK 


— 

Link 

cH  (Packet/Sec) 
i 

1,12 

50 

► V 

2 ,3 

80 

4,9 

70 

5,6 

45 

7,8 

50 

10,11 

70 

Table  7.3 
ROUTING  TABLE 


— 

Source 

Stations 

Destination 

Stations 

Route 

A 

D 

2,4 

A 

E 

2,5 

B 

C 

12,2 

B 

E 

11,8 

C 

B 

4,10 

D 

A 

9,3 

E 

A 

6,3 

E 

B 

7,10 

Table  7.4 

EXTERNAL  ARRIVAL  RATE  (PER  SECOND) 


Node 

Message 

Arrival 

Rate 

Packet 

Arrival 

Rate 

A 

12 

60 

B 

16 

80 

C 

16 

80 

D 

16 

80 

E 

16 

80 

L ji 


I 


I 


Table  7.5a 

MEAN  QUEUE  LENGTH  UNDER  MESSAGE  SWITCHING 


Server 


Simulation 


Method  P 


Method 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 


0.416* 

3.895  ± 0.080 
2.732  ± 0.092 
3 .356  ± 0.081 
3 .225  ± 0.076 
7.289* 

3 .680* 

3 .606  ± 0.165 
1.257* 

5.488  + 0.225 
1.257* 
3.680* 


0.416 
3.947 
2.733 
3 .367 
3.210 
7.289 
3 .680 
3 .654 
1.257 
5.392 
1 .257 
3 .680 


0.536 
4 .272 
3 .033 

3 .681 
3.521 
7.644 

4 .000 

3 .974 
1 .486 
5.73  5 
1 .486 

4 .000 


B 


Method  A 


0.417 


3.862 
2.646 
3 .300 
3.141 
7.210 
3 .617 
3 .537 
1.230 
5 .265 
1.230 
3 .617 


* 

Exact . 


Table  7.5b 

SQUARED  COEFFICIENTS  OF  VARIATION  OF  INTERDE- 
PARTURE TIME  UNDER  MESSAGE  SWITOIING 


Server 

Method  P* 

Method  B* 

Method  A 

1 

0.982 

0.982 

0.8 

2 

0.868 

0.864 

0.8 

3 

0.888 

0.878 

0.8 

4 

0.877 

0.875 

0.8 

5 

0.879 

0 . 876 

0.8 

6 

0 .842 

0 .842 

0.8 

7 

0.872 

0.872 

0.8 

8 

0.872 

0.869 

0.8 

9 

0.935 

0.935 

0.8 

10 

0.853 

0.848 

0.8 

11 

0.935 

0.935 

0.8 

12 

0.872 

0.872 

0.8 

r 



f 

1 

Table  7.6a 

MEAN  QUEUE  LENGTH  UNDER  PACKET  SWITCHING 

1 

Server 

Simulation 

1 

0.364* 

2 

2.399  ± 0.094 

3 

1.621  + 0.051 

4 

2.165  ± 0.061 

5 

2 .027  + 0.063 

6 

4.444* 

7 

2.400* 

8 

2.301  i 0.112 

9 

0.952* 

10 

2.881  ± 0.116 

11 

0.952* 

12 

2.400* 

Method  P Method  B 


2.962 
0.952 
2 .400 


Method  A 


Table  7.6b 

SQUARED  COEFFICIENTS  OF  VARIATION  OF  IN’TERDE 
PARTURE  TIME  UNDER  PACKET  SWITOTING 


mm 


servers  or  the  external  Interarrival  time  distributions  not  only  have 
large  coefficients  of  variation  but  also  are  type  B hyperexponential 
distributions,  the  diffusion  approximation  will  not  work.  This  is  sim- 
ply l^cause  the  arrival  processes  of  the  servers  receiving  jobs  from 
those  servers  with  high  coefficients  of  variation  of  service  times  will 
have  high  coefficients  of  variation.  However,  even  if  the  service  time 
distributions  with  high  coefficients  of  variation  are  type  A hyperexpo- 
nential distributions,  the  diffusion  approximations  under  decomposition 
techniques  still  may  not  be  satisfactory  for  certain  network  topology. 

The  most  noteworthy  networks  of  this  type  are  network  with  feedback 
loops,  especially  self  loops.  Let  us  take  a second  look  at  the  network 
in  Fig.  7.1.  When  the  coefficients  of  variation  of  the  service  times 
in  both  servers  are  not  large,  the  decomposablllty  of  the  network  seems 
to  be  acceptable  from  the  results  obtained  by  various  diffusion  approx- 
imations in  [6]  and  [12]. 

Now,  let  the  service  time  distribution  of  the  first  server  be  hy- 
perexponential distribution  with  the  following  parameter; 

w = 0.029249 
= 30.000187 
1^2  = 0.0232795 

Recall  the  hyperexponential  density  function  has  the  form  (w/M  ) e + 

"X/M  ^ 

((1-w)/M2)  e 2.  This  hyperexponential  distribution  has  mean  0.9  and 

sq.  coeff.  of  variation  64.  One  of  its  branches  has  mean  very  close  to 

zero.  Let  the  second  server  have  exponential  service  time  with  mean 

0.84.  Furthermore,  let 


0^  = 0.5 
02  = 0.5 
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The  simulation  result  of  the  mean  queue  length  in  the  first  server 
with  95^  confidence  interval  is  165+15.  The  mean  queue  lengths  obtained 
by  diffusion  approximation  are  around  294  under  methods  A and  13  and  264 
under  method  P.  That  is  to  say,  all  methods  tend  to  overestimate  the 
mean  queue  length.  We  now  try  to  give  a reasonable  explanation  to  this 
anomaly.  Recall  the  arrival  rate  to  the  first  server  can  be  calculated 
by  solving  the  system  of  linear  equations  (7.3).  After  simple  manipula- 
tion, we  get  = 1 . Notice  the  external  arrival  rate  is  0.5,  so 

one-half  of  the  arrivals  to  the  first  server  is  from  the  feedback  path 
through  the  second  server.  When  the  first  server  encounters  a long  ser- 
vice time,  the  arrival  process  from  the  second  server  will  be  shut  down 
after  the  second  queue  becomes  empty.  That  is  to  say,  the  arrival  rate 
is  effectively  0.5  instead  of  1 during  the  later  period  of  the  long  sci’- 
vice  time.  However,  in  the  diffusion  approximations,  the  arrival  rate 
to  the  first  server  is  always  1.  This  is  the  reason  why  the  mean  queue 
length  is  overestimated  under  diffusion  approximation.  Tlie  corrolations 
of  the  seri’ice  stations  become  very  serious  as  the  coefficients  of  vari- 
ation of  service  time  distributions  become  large  in  this  typo  of  network, 
and  using  the  ordinary  way  to  estimate  the  diffusion  parameters  is  not 
sufficient  to  account  for  this  sort  of  correlations.  Wliat  will  happen 
if  the  service  time  distribution  is  a type  B hyperexponential  distribu- 
tion function  indicated  in  Section  5?  A simulation  has  been  conducted 
when  the  first  server  has  a two  stage  hyperexponential  distribution  with 
the  following  parameters: 


w = 0.000126 
= 500 

= 0.843486 


i ■ 


I 


The  second  server  still  has  exponential  distribution.  The  network  topol- 
ogy and  the  traffic  intensity  is  the  same  as  the  previous  example.  Hence, 
the  moan  queue  length  predicted  by  diffusion  approximation  is  still  the 
same.  But  the  951^  interval  estimation  obtained  liy  simulation  is  9(3  ->  27, 
which  is  quite  small,  as  expected.  The  broad  wicltli  of  the  confidence 
interval  is  due  to  the  heavy  traffic  condition  and  the  closeness  of 
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to  0.  The  number  of  arrivals  in  the  simulations  is  around  10®,  hence 
the  point  estimation  should  be  acceptable. 

So,  the  idea  of  decomposing  a network  of  queues  into  separate  single 
server  queueing  systems  does  not  seem  to  be  always  feasible.  In  certain 
network  topologies,  such  as  network  with  feedback  loops,  especially  self 
loops,  there  is  a dependence  of  the  arrival  process  of  each  service  cen- 
ter in  the  feedback  path  on  its  departure  process.  If  the  service  center 
has  a self  loop  which  contributes  a large  portion  of  arrivals,  the  effect 
of  this  dependency  becomes  very  serious  as  the  coefficient  of  variation 
of  the  service  time  deviates  significantly  from  1.  Using  any  of  the  three 
methods  cited  above  to  estimate  the  coefficient  of  variation  of  the  in- 
terarrival time,  the  decomposition  technique  can  not  reflect  this  depen- 
dency into  the  estimated  parameter. 

To  be  more  specific,  let  us  consider  method  P.  This  fact  can  be 
observed  from  Fig.  7.2b  where  the  two  service  centers  with  rate  are 

actually  the  same  but  are  represented  as  two  different  service  centers. 
Clearly,  the  dependence  among  the  two  is  not  reflected  in  the  estimated 
parameters.  Nevertheless,  the  self  loop  problem  can  be  solved  by  treat- 
ing the  server  and  its  self  loop  as  a single  entity  as  we  did  in  Section 
6.  That  is  to  say,  we  first  eliminate  all  self  loops  in  the  network  by 
replacing  each  server  with  a self  loop  by  an  equivalent  server  witliout  a 
self  loop  as  In  Fig.  7.4  and  then  apply  the  decomposition  technique  if 
possible.  Ijet  P'  be  the  contribution  to  the  diffusion  parameter  p from 
the  server  and  its  self  loop  if  any.  In  Table  7.7,  we  tabulate  the  ap- 
proximate values  of  p'  under  methods  P*  and  A*  when  the  self  loop  is 
not  eliminated  and  the  correct  value  under  the  equivalent  server  without 
a self  loop.  The  error  terms  under  the  two  approximation  methods  are 
proportional  to  0^(0“^  +0).  This  explains  why  direct  decomjx>s  ition  docs 
not  work  for  a strong  self  loop  under  a large  coefficient  of  variation 
of  service  time  even  if  the  distribution  is  type  A.  Besides  the  anomal- 
ies due  to  type  B distributions,  the  problem  still  not  solved  is  strong 
feedback  loops  which  are  not  self  loops  under  large  coefficients  of  vari- 
ation of  service  times.  Although  the  analysis  is  greatly  simplified  when 
decomposition  does  work,  decomposition  is  not  a panacea.  We  should  bo 
careful  about  the  decomposability  of  the  queueing  network  and  all  tlic 
distributions  involved. 
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Server  with  a self  loop 
Variance  of  service  time 


Equivalent  server  without  a self  loop 
Variance  of  service  time;  c^/Cl-©)  + 

e/(u(i-e))^ 


Fig.  7.4.  SERVER  WITH  A SELF  LOOP  AND 
ITS  EQUIVALENT  REPRESENTATION  WITHOUT 
A SELF  LOOP. 


Table  7.7 


UNDER  VARIOUS  METHODS 


With  a Self  Loop 


Without  a 
Self  Loop 


Motluxl  P 


Metliod  A 


8.  THE  SERVICE  CENTER  WITH  A QUEUE  DEPENDENT 
SERVICE  RATE  OR  ARRIVAL  RATE 


In  this  section,  we  consider  the  case  where  a service  center  has 
queue  dependent  service  rate.  The  conventional  G/M/m  queueing 
system  is  a special  case  of  this  class  of  service  centers.  The  service 
rate  of  the  m-server  queueing  system  can  be  expressed  as 


for 

i < m 

^‘i  = 

* mp 

for 

i > m 

where  i is  the  number  of  customers  in  the  queue . 

In  computer  system  modeling,  a more  general  than  the  conven- 

tional one  cited  above  is  often  needed.  Consider  the  performance  of  a 
tightly  coupled  computer  system.  The  total  service  rate  of  CPU’s  does 
not  increase  as  a linear  function  of  the  number  of  CPU's  in  the  system 
due  to  the  memory  interference  among  different  CPU's. 

When  the  service  rate  is  queue  dependent,  the  diffusion  parameters 
also  become  queue  dependent  and  the  diffusion  equation  becomes  harder  to 
solve.  We  first  need  to  determine  the  values  of  the  diffusion  parameters 
at  those  integer  points  and  then  propose  a reasonable  way  to  interpolate 
their  values  in  between  Integers.  The  infinite  capacity  case  is  first 
considered.  We  further  assume  that  will  keep  constant  for  i ' m, 

as  in  the  conventional  multiserver  case  with  m servers.  Similarly,  for 
the  arrival  rate 

The  values  of  the  diffusion  parameters  at  those  integer  points  are 
defined  as 

(8.1) 

(«!  = C^X^  + C^gp. 

and  as  usual,  is  defined  to  be 


r . 
1 
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where  C and  C are  the  squared  coefficients  of  variation  of  the 

cl  S 

interarrival  time  and  service  time  distributions,  respectively.  When 

the  service  rate  is  fixed,  we  set  g equal  to  the  traffic  intensity  of 

the  queueing  system.  But,  for  a server  with  queue  dependent  service 

rate  or  arrival  rate,  the  appropriate  value  of  g is  not  very  clear. 

From  experimental  results,  it  seems  to  be  that,  if  we  set  g equal  to 

1 for  the  case  where  C < 1/2  and  C <1  and  A /u  otherwise , the 

a — s — mm 

approximation  will  be  more  accurate  in  general. 

There  are  at  least  two  different  ways  to  interpolate  the  value  of 
Q(X)  and  P(X)  in  between  integers. 


Method  1.  Interpolation  by  Step  Functions  (see  Fig.  8.1) 


Q(X)=<'a,  k-|<X<k+|  and  2 < k ^ m 


(8.2) 


la 


X > m - - 


Similarly  for  3(X) . 


Method  2.  Interpolation  by  Linear  Functions  (see  Fig.  8.2) 


X _ 1 


/(X)  = + (a  -rO(X-k)  k<X<k+l  and  1 k m-1 

'Itk+lk  — — - 


X m 


(8.3) 


or  lot  a(X)  = for  X < 1 to  make  the  impulse  term  6(X-1) 

in  the  diffusion  equation  easier  to  handle.  Similarly  for  . 


Again  , 
cf)nd  i t ion  . 


wo  use  Feller's  elementary  return  process  to  handle  boundary 
The  diffusion  equation  satisfied  by  the  pj’obability  density 
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Fig.  8.1.  interpolation  OF  a(X)  USING  STEP 
FUNCTIONS . 


Fig.  8.2.  INTERPOLATION  OF  0;(X)  USING  LINEAR 
FUNCTIONS. 


V'* 


► 


function  P(X)  of  the  approximate  queue  length  process  has  the  follow- 
ing form  under  steady  state 


i a(X)  P(X)  - ^ p(X)  P(X)  = -AqM^6(X  - 1)  (8.1) 

d X 


with  boundary  conditions 


and 


lim 
X->  0 


^ a(X)  P(X)  - p(X)  P(X)  = 

2 dX  u 1 


P(0)  = 0 


where  is  the  probability  that  the  queueing  system  is  idle. 

Under  the  second  interpolation  method,  numerical  integration  is  re- 
quired to  estimate  the  queue  length  distribution.  For  the  conventional 
multiserver,  the  broken  line  in  Fig.  8.2  becomes  a straight  line  and  the 
complexity  of  the  problem  is  simplified,  but  numerical  integration  is 
still  needed.  Hence,  for  better  mathematical  tractabil ity , we  adopt  the 
first  interpolation  method.  Halachmi  and  Franta  [26,25]  have  applied 
diffusion  approximation  to  conventional  multiserver  queueing  systems  us- 
ing the  second  interpolation  method  and  reflection  boundary  for  both  in- 
finite and  finite  population  models,  respectively.  The  results  from  both 
methods  seem  to  be  quite  close  in  the  few  cases  examined. 

After  solving  the  differential  equation,  we  got 


^ V / 


(X-l)r, 


M^d  e 


P(X) 


M^d  o 


M.^d  e 


8 ,+(X-m+-)r 

m-1  2 m 


for  0 X 1 


for  1 < X < - 


(8.5) 


for  k - i X h + i and  2 _ k _ m - 1 


1 

for  X ■ m - - 
81 


J 
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j: 

i 


I 

I 


■ 


J 


where 


k 

V "I 

'k  ^k  - T 


s.  = 


for  1 < k < m - 1 


(«.6) 


and 


i 

F - 


1«‘ 


The  unknown  constant  can  be  determined  by  the  fact  that  total 

probability  must  sum  up  to  1,  i.e., 


t ^ » 

I 


f 

•'r\ 


P(X)  dX  + = 1 


After  simplification,  we  get 


i 


\ 


In  Table  8.1,  we  compare  the  mean  queue  lengths  obtained  under  dif- 
fusion approximation  with  analytic  results  for  the  M/M/m  system  with 
m = 2 , 3 , 1 , 5,6,  7 , and  8,  when  p ^ A/mu  = 0.95  and  0.85,  respectiv'ely . 
The  approximation  is  very  accurate.  Then,  we  compare  the  conditional 
moan  queue  length  of  the  external  queue  (given  that  external  queue  ex- 
ists, i .e . , number  of  iobs  in  the  system  is  larger  than  m)  with  the  an- 
alytic result  for  the  G/M/m  queueing  system.  Both  the  analytic  result 
and  diffusion  approximation  on  the  conditional  mean  external  queue 
length  of  the  G/M/m  queueing  system  are  independent  of  m . After  simple 
manipulation,  we  can  got  the  conditional  mean  external  queue  length  un- 
der  diffusion  appro,ximat ion  which  is  1/(1  -e  ).  The  exact  result  is 
1^(1  -a)  whore  a is  defined  in  Section  3.  In  Table  8.2,  we  compare 
the  conditional  mean  external  queue  length  obtained  under  diffusion  ap- 
proximation with  the  analytic  result  for  E^/M/m  and  E^/M/m  c|ueuoing  sys- 
tems when  p = 0.95,  0.90,  0.85,  0.80,  and  0.75. 

When  the  arrival  process  has  hypcrexponontia 1 distribution,  again, 
the  conditional  mean  external  queue  length  can  vary  over  a wide  range. 
Nevertheless,  for  type  A hyperexponential  distribution,  diffusion  ap- 
proximation can  still  be  applied  as  before.  Tables  8.3  and  8.1  tabulate 
the  diffusion  approximations  and  analytic  results  under  different  values 

of  M for  C = 2,  1,  8,  16,  and  32,  when  p = 0.95  and  0.85,  respec- 

2 a 

tively.  In  both  cases,  A is  equal  to  1. 

In  Table  8.6,  we  consider  the  case  whore  is  an  arbitrary  func- 

tion of  i and  A is  constant.  Tlie  result  is  again  satisfactory, 
i 

We  now  consider  the  closed  two  server  queueing  network  in  Fig.  8.3a 
which  can  be  Interprctefl  as  the  CPlf/DTU  model,  as  noted  earlier. 
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Table  8.1 

MEAN  QUEUE  LENGTHS  FOR  M/M/m  SYSTEM  WHEN  p = 0.85  AND  0.95 
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Generalization  to  the  model  in  Fig.  8.3b  follows  the  same  idea  cited  in 
Section  6.  That  is  to  say,  we  first  compute  the  mean  and  variance  of 
the  service  time  of  the  equivalent  server  without  a self  loop.  The  only 
difference  from  before  is  that  the  quantities  now  are  queue  dependent. 
After  replacing  each  server  and  its  self  loop  by  an  equivalent  server 
without  a self  loop,  the  model  in  Fig.  8.3b  reduces  to  that  in  Fig.  8.3a. 
The  number  of  jobs  in  the  system  is  assumed  to  be  K.  The  queue  depen- 
dent service  rate  of  CPU  and  DTU  are  denoted  by  p and  A \dien 

X K"  X 

there  are  i jobs  in  the  CPU  queue.  The  diffusion  parameters  are  de- 
fined to  be 


h ■ \-i  - "i 

“i  = ♦ '=,“1  1 < * £ K 

and , as  before , 


(S.9) 


Notice  the  definition  is  similar  to  the  previous  one  with  g = 1.  Again, 
there  are  at  least  two  different  ways  to  interpolate  the  value  of  a(X) 
and  p(X)  in  between  integers.  The  only  difference  is  that  now  we  have 
two  boundaries.  We  still  adopt  the  interpolation  method  using  step 
functions  for  simplicity.  To  be  more  precise,  we  define 


0<X<- 


a(x)  = < O', 


i -|<X<1  +1 


2 <KK  -2 


(8.10) 


.Vl 


K - 2 < X < K 


Similarly  for  3(X) . 


The  diffusion  equation  now  has  the  following  form 


I a(x)  p(x)  - ^ p(x)  p(x)  = -AqM^6(x  - 1)  - m^m^SCx  - m + i) 
dX 


lim  i ^ a(X)  P(X) 
X^  0 2 dX 


P(X)  P(X)  = 


(8.11) 


lim  I ^ a(X)  P(X) 
X^K  2 


- p(x)  P(X)  = 


where  Mg  is  the  probability  that  the  CPU  queue  has  K jobs . 
After  solving  the  diffusion  equation  (8.11),  we  get 


Vi  ' 


for  0 < X < 1 


(X-l)r, 


M^d  e 


for  1 < X < - 


P(X)  = 


Mj^d  e 


Mj^d  e 


(8.12) 


for  1 - i < X < max  (K  - 1,  i + ^) 


\-2 

\ - 1 


and  2 < i < m - 1 


for  X > K - 1 


where 


“l^^K-1  \-2‘*’2^K-l 


(8.13) 


and  is  defined  as  before  in  (8.6). 

The  unknown  constant  can  be  determined  by  the  fact  that  total 

probability  must  sum  up  to  one,  i.e., 
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1 


P(X)  dX  + M,  + M = 
X ^ 


After  simplification,  we  get 


K-2 


,5 


(8.14) 


+ d ^(r^^_i-l)  e ^ ^+1^  ^^K-1 


1-1 


/ 


Let  be  the  probability  that  i jobs  are  in  the  CPU.  We  define 


rto  =Mi 


■ /„ 

C-.  - r 


3/2 


P(X)  dX 
i+1/2 


P(X)  dX 


1/2 


IC.  1 = I P(X)  dX 

\-3/2 


\ \ = 


= M„ 


In  Table  8.6,  we  compare  the  mean  queue  length  at  the  CPU  when  both 
CPU  and  DTU  have  exponential  service  time  distributions,  and  furthermore 
the  CPU  is  modeled  as  a traditional  m-server  under  fixed  degree  of  multi- 
programming K.  In  Table  8.7,  the  case  where  service  rate  of  CPU  is  an 
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Table  8.6 

MEAN  QUEUE  LENGTH  WHEN  CPU  IS 
MODELED  AS  CONVENTIONAL  m SERVER 


y ' 

m 

— 

K 

p = 0. 

95 

p = 0. 

85 

Diffusion 

Exact 

Diffusion 

Exact 

» 4 

2 

4 

2.28 

2.13 

2.11 

1.94 

F 

3 

2.00 

1.92 

1.91 

1.81 

i ' 

3 

5 

3.09 

2.89 

2.88 

2.64 

r : 

1 

4 

4 

2.81 

2.69 

2.70 

2.55 

4 

7 

4.38 

4.14 

4.02 

3.73 

L ' 

5 

5 

3.64 

3.49 

3.48 

3.31 

5 

8 

5.19 

4.94 

4.78 

4.48 

6 

6 

4.47 

4.31 

4.27 

4 .08 

. 

6 

10 

6.48 

6.22 

5.89 

5.56 

r 

7 

7 

5.30 

5.14 

5.07 

4.87 

r 

h 

7 

11 

7.31 

7.05 

6.66 

6.33 

8 

8 

6.15 

5.98 

5.87 

5.66 

8 

^ - — 

12 



8.15 

7.89 

7.44 

7.11 

I 
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I arbitrary  function  of  queue  length  and  that  of  DTU  is  constant  is  con- 

1 • 

sidered.  Both  the  mean  queue  length  and  utilization  at  the  CPU  is 
tabulated.  Again,  the  result  is  quite  satisfactory. 


9.  C»NCLUSION 


Diffusion  approximation  is  an  attractive  means  of  approximating  tlio 
performance  of  queueing  systems.  In  this  paper,  we  not  only  assess  the 
accuracy  of  diffusion  approximation  but  also  its  limitation  and  applica- 
ble range.  Modern  computer  systems  are  so  complicated  that  oversimpli- 
fied models  may  not  predict  any  useful  results.  Realistic  models  often 
are  not  analytically  tractable.  Finding  approximate  solutions  or  upper 
bounds  and  lower  bounds  of  the  solutions  is  the  only  means  to  handle  more 
complicated  problems  short  of  simulation.  Under  heavy  traffic  conditions, 
simulation  converges  very  slowly  and  diffusion  approximation  seems  to  bo 
the  most  attractive  way  to  solve  the  problem.  Nevertheless,  diffusion 
approximation  is  not  a panacea,  it  does  have  a limitation.  This  limita- 
tion has  been  overlooked  in  the  past.  Substantial  effort  has  been  de- 
voted in  this  paper  to  identify  the  conditions  where  diffusion  approxi- 
mation can  obtain  accurate  estimates.  We  must  be  careful  with  these 
conditions  when  applying  diffusion  approximation. 

In  Table  9.1,  we  classify  the  single  seihrer  queueing  systems  accord- 
ing to  their  coefficients  of  variation  of  service  times  and  interarrival 
times,  and  point  out  the  diffusion  approximation  technique  which  seems  to 
be  most  accurate  according  to  our  analysis.  The  superiority  of  method 
the  proposed  method,  should  be  very  apparent.  When  is  larger  than 

one,  the  mean  queue  length  may  vary  over  a wide  range  even  if  the  first 
two  moments  of  Interarrival  time  are  kept  constant.  Diffusion  approxima- 
tion is  applicable  under  the  condition  that  the  high  variation  of  inter- 
arrival  time  is  due  to  a great  number  of  short  Interarrival  times  instead 
of  a few  very  long  interarrival  times.  Case  studies  have  been  conducted 
on  2-stage  hyperexponential  distributions  wblch  are  widely  used  in  com- 
puter system  modelling.  A similar  anomaly  is  observed  in  two  server 
closed  queueing  networks,  often  referred  to  as  .CPU/DTU  models,  when  the 
service  time  of  any  server  has  a large  coefficient  of  variation.  Again, 
a similar  regularity  condition  on  service  time  distributions  is  required 
in  order  for  the  diffusion  approximation  to  be  applicable.  Although 
method  B does  not  yield  the  best  performance  \»^en  applying  it  to  approx- 
imate the  single  server  system,  it  is  Indeed  a nice  way  to  approximate 
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Table  9.1 

RECOMMENDED  DIFFUSION  APPROXIMATION  METHOD  FOR  SINGLE  SERVER  SYSTEM 

i . 


Cg  Close  to  or 
Less  than  0.5 

C Close  to  1 
s 

C > 1 
s 

C < 1 
a 

Method  A 

Method  P 

Method  P 

C ^ 1 
a 

Method  P 

Method  P 

Method  P 

C > 1* 
a 

1 

Any  Method 

Method  P 
(?) 

The  high  coefficient  of  variation  of  interarrival  time  must 
be  due  to  a large  number  of  short  Interarrival  times  in- 
stead of  a few  very  long  interarrival  times. 


the  two  server  closed  queueing  network.  When  the  coefficient  of  varia- 
tion of  the  CPU  service  time  is  small,  method  G1  has  similar  performance. 
As  the  coefficient  of  variation  of  the  CPU  service  time  increases  , method 
P becomes  somewhat  better. 

For  general  queueing  networks,  an  efficient  way  of  taking  into  ac- 
count the  effect  of  idle  periods  to  estimate  the  coefficient  of  variation 
of  the  arrival  process  at  each  server  when  the  network  can  be  decomposed 
into  separate  single  servers  is  proposed.  For  certain  network  topolo- 
gies, the  arrival  processes  of  some  service  centers  strongly  depend  upon 
their  o\m  departure  processes.  Networks  of  this  type  are  networks  with 
strong  feedback  loops,  especially  self  loops.  When  the  coefficients  of 
variation  of  the  service  times  at  the  service  centers  have  a large  devi- 
ation from  one.  this  sort  cif  queueing  network  can  not  be  decomposed  into 
separate  single  servers  directly.  This  fact  has  been  neglected  in  the 
past.  Nevertheless,  the  self  loop  problem  can  be  solved  by  replacing 
e.ach  seiner  with  a self  loop  by  an  equivalent  server  without  a self  loop, 
'.ftei  eliminating  all  the  self  loops,  we  can  reconsider  the  decomposition 
of  a network.  The  problem  still  not  solved  seems  to  be  networks  with 
-tronp  feertback  loops  v’hich  are  not  self  loops  when  the  coefficients  of 
variation  of  some  service  times  are  large.  Surely,  the  regularity  con- 
diti'ju  th  it  a laige  coefficient  of  variation  of  external  interarrival 
time  or  service  time  of  each  intermediate  server  is  due  to  a lot  of  short 
intcrari’i\;il  times  or  service  times,  respectively,  must  always  hold  in 
order  for  diffusion  approximation  to  be  applicable. 

Finally,  we  consider  the  service  center  with  queue  dependent  service 
rate  or  arrival  rate.  General  queue  dependent  service  rate  is  often  en- 
countercfl  in  computer  system  modeling.  Generalization  to  closed  two  ser- 
v'er  queueing  netwoi-k  where  each  serv'er  may  have  a self  loop  is  also  con- 
s idored . 
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